DOE/ER/03077-178 
MF  -  98 


Courant  Institute  of 
Mathematical  Sciences 

Magneto-Fluid  Dynamics  Division 


Magnetic  and  Drift  Surfaces  in 
Toroidal  Plasma  Equilibria 

Michael  Marcal 


Prepared  under  Contract    DE-AC02-76ER03077 
with  the  U.  S.  Department  of  Energy 

Research  and  Development  Report 

Plasma  Physics 
October  1982 


New  York  University 


UNCLASSIFIED 

New  York  University 

Courant  Institute  of  Mathematical  Sciences 

Magneto-Fluid  Dynamics  Division 


MF-98  DOE/ER/03077-178 

MAGNETIC  AND  DRIFT  SURFACES  IN 
TOROIDAL  PLASMA  EQUILIBRIA 

Michael  Marcal 
October  1982 


U.    S.    Department   of  Energy 
Contract   No.    DE-AC02-76ER03077 


UNCLASSIFIED 


\ 


-111- 
CONTENTS 

Abstract  v 

Acknowledgements  vi 
Chapter  1 

A  general  discussion  of  the  problem                               1 

Chapter  2 

The  computational    scheme   of    Betancourt  and  Oarabedian  12 

Chapter    3 

A   methDd    to   study   magnetic    surfaces    in    toroidal    plasmas 

by   following    field   lines  20 

(a)  The    representation   of    the    nagnetic    field 

(b)  The   field   line  equations 

Chapter    4 

The  application   of   the    method    to   study   drift  surfaces                                  29 

(a)  The    guiding  center   equations  29 

(b)  Scaling    the  orbit   eqiiations  37 

Chapter    5 

A   discussion   of   some    examples   and    results  ^^ 

(a)  The    equilibria    studied  '^O 

(b)  Examples  "^2 
Example  1.  M^ignetic  surfaces  42 
Example  2.      The   effect    of    not    preserving    the  divergence-free 

condition    on    solenoldal    vectors  5L 

Example   3.      Drift   surfaces    and    tVie    dependence    on    the    gyroradius  56 


-IV- 

Exaraple  4.      Trapped   particle   orbits  63 

(c)   Discussion  68 

Chapter   6 

The  computer    programs   COEFF  and  LINOIB  70 

(a)  The   operation  of   the   programs  70 

(b)  The  FORTRAN  listings   of    COEFF  and  LINORB  75 

References  128 


ABSTRACT 

We  have  developed  a  method  for  following  magnetic  field  lines  and 
guiding  center  orbits  of  isolated  charged  particles  in  toroidal  plasma 
configurations  with  no  spatial  symmetries.  The  heart  of  the  method 
lies  in  a  divergence-free  approximation  to  a  numerical  equilibrium 
magnetic  field  B_  that  has  been  computed  over  a  grid  of  points  by  the 
computer  code  of  Betancourt  and  Garabedian  [3].  Considerations  from 
the  KAM  theory  for  Hamiltonian  systems  and  the  results  of  our  studies 
lead   us    to  believe    that    the   solenoidal    property    of  _B, 

V    .    B   =   0 
must      be      preserved      to      be     able      to     follow  trajectories    in    the    torus 
adequately . 

The  Be tancourt-Garabedian  code  is  first  used  to  compute  an 
equilibrium  over  a  rectangular  grid  in  a  unit  cube  0  <^  s  _<  1, 
0<u<l,  0<v<l  into  which  the  plasma  or  vacu\jm  region  has  been 
assumed  to  be  mapped.  In  this  coordinate  system,  s  denotes  a 
radial-like  flux  variable,  and  u  and  v  are  angular  variables  with  unit 
periods.  Implicit  in  their  formulation  is  the  hypothesis  of  the 
existence  of   a   nested    family   of    flux   surfaces. 

Our  approximation  to  the  numerical  magnetic  field  output  by  their 
code    is   based    on   a    representation    of    B  as 

Vii)  X  Vu  +  Vx  X  Vv 
where  the  functions  ij;  and  x  ^re  series  of  polynomials  in  s,  e^  '^  and 
gi-Tiv  ^l^\■^  coefficients  to  be  determined  by  a  simple  least  squares  fit 
to  the  discrete  magnetic  field,  Oiir  representation  does  not  assume 
t>iat  the  flux  surfaces  are  nested,  as  examples  will  exhibit.  Moreover, 
the         representation         leads         to        very      simple      systems      of      ordinary 
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differential  equations  that  describe  the  field  lines  and  guiding  center 
orbits.  We  have  employed  Boozer's  formulation  of  the  guiding  center 
equations    [9], 

The  method  has  been  applied  to  study  magnetic  flux  surfaces,  drift 
surfaces,  and  also  the  orbits  of  trapped  particles  in  finite  beta 
plasma  configurations  with  stellarator  geometry.  We  have  included  the 
FORTRAN  listings  of  two  computer  programs  that  implement  our  method. 
These  have  been  written  to  be  fully  compatible  with  the  equilibrium  and 
stability  code   of    Betancourt  and  Garabedian. 
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1.    General   Discussion  of   the  Problem 

The  ability  to  maintain  a  plasma  in  equilibrium  is  a  central  goal 
of  magnetic  fusion  energy  research.  The  plasma  has  to  be  confined  by  a 
strong  magnetic  field  in  a  toroidal  vessel  long  enough  for  the 
controlled  release  of  energy  to  take  place.  In  particular,  the 
equilibrium  problem  for  closed  systems  is  to  specify  a  magnetic  field 
in  toroidal  geometry,  consistent  with  Maxwell's  equations,  so  as  to 
ensure  the  long  term  confinement  of  charged  particles  within  the 
vessel. 

The  containment  of  a  plasma  in  such  geometry  rests  on  the 
assumption  of  the  existence  of  a  nested  set  of  toroidal  magnetic  flux 
surfaces  lying  within  the  confining  vessel.  \  consequence  of  the 
equations  of  ideal  magnetohydrodynamics ,  which  describe  the  plasma  as  a 
continuous  fluid,  is  that  plasma  cannot  flow  across  flux  surfaces  [16]. 
Therefore,  the  presence  of  a  nested  family  of  such  surfaces  guarantees 
the   confinement    of    the    plasma   as    a   bulk.. 

The  origin  of  a  flux  surface  is  easily  visualized  by  following  a 
magnetic  field  line  trajectory  many  times  around  the  torus.  If  the 
field  line  remains  within  the  torus,  there  are  three  possibilities  for 
its  behavior.  The  first  is  that  the  field  line  may  close  on  itself 
after  a  finite  number  of  transits  of  the  torus.  The  second  is  that  the 
field  line  may  lie  on  a  closed  toroidal  surface  which  it  covers 
ergodically.  This  means  that  the  field  line,  sooner  or  later,  comes 
arbitrarily  close  to  any  point  of  that  surface;  such  a  surface  is 
called  a  magnetic  surface.  The  flux  of  the  magnetic  field  across  a 
plane  bounded  by  a  magnetic  surface  and  transverse  to  it  is,  by  an 
application      of       the      divergence    theorem,    a    constant    independent    of    the 
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transverse  plane  and  depends  only  on  the  magnetic  surface.  Hence,  a 
magnetic  surface  is  also  called  a  flux  surface.  The  third  possibility 
is  that  the  field  line  trajectory  ergodically  fills  a  region  of 
positive  volume  within  the  torus.  This  is  probably  the  most  general 
situation  for  an  arbitrary  magnetic  field  whose  lines  of  force  remain 
within  the  torus:  the  occurrence  of  magnetic  surfaces  is  a  rather 
special  feature  of  certain  magnetic  fields.  In  fact,  the  existence  of 
such  surfaces  has  only  been  proved  rigorously  when  the  magnetic  field 
is  assumed   to   possess   a  high   degree   of   symmetry    [2,16]. 

The  nestedness  requirement  of  the  flux  surfaces  is  also  believed 
to  be  necessary  for  the  longtime  confinement  of  charged  particles  in 
the  plasma.  An  analysis  of  the  ordinary  differential  equations  that 
describe  the  orbit  of  a  charged  particle  in  a  magnetic  field  that  has  a 
strong  toroidal  component  reveals  that  the  field  line  trajectory  is  the 
lowest  order  approximation  to  the  orbit.  Hence,  the  nested  character 
of  the  toroidal  flux  surfaces  would  imply  particle  confinement,  at 
least    to  lowest    order. 

The  next  higher  order  approximation  to  the  particle's  orbit  is 
called  the  trajectory  of  the  guiding  center.  The  actual  motion  of  the 
particle  consists  of  a  rapid  gyration  about  an  instantaneous  center, 
which  generally  follows  the  path  of  a  field  line  but  which  may  also 
drift  across  field  lines  after  a  long  enough  time.  This  drift  motion 
of  the  guiding  center  occurs  because  of  the  curvature  of  the  field  line 
and   the  variation  of   the   magnetic    field   strength  along    the   orbit. 

If  a  particle  is  able  to  traverse  the  torus  many  times,  the  orbit 
of  its  guiding  center,  being  a  perturbation  of  a  field  line,  should 
sweep      out     a   toroidal   drift    surface    which   would  be   close   to  a    magnetic 
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surface,  if  the  latter  exists.  It  seems  reasonable  to  expect  good 
confinement  of  the  plasma  when  there  are  nested  flux  surfaces  and 
nearby    nested   drift    surfaces. 

This  last  condition  is  by  no  means  necessary  or  sufficient  for 
confinement.  First  of  all,  the  behavior  of  a  single  particle  not 
subject  to  collisions  can  surely  not  be  an  adequate  description  of  the 
collective  behavior  of  particles  in  a  plasma.  Second,  the  orbit  of  the 
guiding  center  is  an  approximate  description  of  the  particle's  true 
path:  including  higher  order  terms  may  reveal  radically  different 
behavior  than  that  predicted  on  the  basis  of  guiding  center 
calculations.  Related  to  this  is  the  fact  that  small  perturbations  of 
the  magnetic  field  may  have  a  drastic  effect  on  the  geometry  of  flux 
surfaces  and  drift  surfaces.  For  example,  a  nested  family  of  toroidal 
magnetic  flux  surfaces  may  be  broken  up  into  regions  where  these 
surfaces  have  quite  complicated  geometries;  or  zones  of  ergodic 
behavior  may  arise  where  the  field  lines  do  not  sweep  out  any  surface. 
Yet  even  in  such  cases,  it  is  not  possible  to  predict  poor  confinement, 
since  such  regions  may  be  relatively  small  and  surrounded  by  closed 
toroidal  surfaces. 

The  object  of  the  present  work  is  to  describe  a  method  of 
representing  the  magnetic  field  in  equilibrium  plasmas  so  that  the 
magnetic  field  lines  and  the  guiding  center  orbits  of  isolated  charged 
particles  may  be  numerically  integrated  to  obtain  the  flux  surfaces  and 
drift  surfaces.  We  shall  be  interested  in  equilibrium  configurations 
of  the  plasma  that  possess  no  spatial  symmetries.  The  numerical 
procedure  and  computer  code  of  Betancourt  and  Harabedian  [3]  is  the 
best   available   numerical    tool    to   search   for  and   analyze      the      stability 
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properties  of  such  equilibria,  and  our  method  is  intimately  linked  to 
this  code.  The  latter  provides  a  numerical  equilibrium  magnetic  field 
defined  in  flux  coordinates.  We  then  approximate  this  field  by  a 
method  that  does  not  assume  the  nestedness  of  the  flux  surfaces  which 
such  a  flux  coordinate  system  entails.  The  study  of  flux  surfaces  and 
drift  surfaces  by  numerical  integration  is  not  new,  but  most  previous 
efforts  have  applied,  strictly  speaking,  only  to  vacuum  configurations 
with  no  plasma  present.  The  significance  of  our  work  is  that  the 
method  is  applicable  to  plasma  and  vacuum  configurations.  Moreover,  in 
our  model  of  the  magnetic  field  _B,  the  solenoidal  character  of  the 
field  is  preserved  exactly.  The  importance  of  the  divergence-free 
condition  on  _B  for  following  field  lines  and  particle  orbits  has  not 
been   emphasized    sufficiently  in   previous   work.      '■■-'  >'   ' 

Gelfand  et  al  [11]  were  probably  the  first  to  look  into  the  nature 
of  flux  surfaces  by  integrating  numerically  the  field  line  equations. 
They  were  led  to  undertake  such  a  procedure  because  the  traditional 
techniques  of  studying  flux  surfaces  by  asymptotic  methods  fail  to 
apply  in  the  vicinity  of  the  magnetic  axis  [8].  Precisely  in  this 
region,  they  found  that  small  perturbations  to  a  magnetic  field  with 
helical  symmetry  in  cylindrical  geometry  would  usually  destroy  a  nest 
of    flux    surfaces. 

Gibson,  in  his  numerical  studies  on  the  nature  of  flux  surfaces 
[12]  and  guiding  center  orbits  [13]  in  £  =  3  stellarator  geometries, 
demonstrated  quite  convincingly  the  effectiveness  of  such  computer 
methods  not  only  as  means  of  checking  theoretical  calculations,  but 
also  as  valid  research  tools  in  their  own  right.  For  the  cases  that  he 
was      Interested      in,    the    magnetic    field  was  assumed    to  be    produced   by  a 
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set  of  filamentary  cur  rent -carrying  coils  surrounding  a  vacuum  chamber. 
It  was  computed  at  each  step  of  the  integration  of  the  field  line  or 
orbit  trajectory  by  means  of  the  Biot-Savart  law  of  magnetostatics. 
This  method  of  using  the  Biot-Savart  law  which  Gibson  exploited  so 
successfully  has  since  become  quite  widespread.  The  use  of  the 
Biot-Savart  law  applies,  strictly  speaking,  only  to  harmonic  (vacuum) 
magnetic  fields.  The  magnetic  field  of  a  plasma  in  equilibrium  is  not 
a  harmonic  vector  field.  So,  any  conclusion  about  the  nature  of  flux 
surfaces  in  equilibrium  plasmas  based  on  calculations  using  vacuum 
fields  with  the  same  geometry  may  be  misleading  even  when  applied  to 
low    pressure   plasmas. 

Another  method  for  tracing  field  lines  in  plasmas  might  involve 
solving  for  the  equilibrium  magnetic  field  on  a  fixed  grid,  and  then 
making  high  order  local  approximations  for  the  magnetic  field  where  it 
is  needed  by  tlie  integration  scheme.  The  disadvantage  of  this 
procedure  is  that  the  solenoidal  character  of  the  magnetic  field  is  not 
preserved  in  the  local  approximation,  or  at  most  only  a  discrete 
version  of  it  is.  It  is  a  contention  of  the  present  paper  that  to 
follow  field  lines  and  particle  orbits  adequately  over  a  long  period  of 
time,  the  solenoidal  character  of  the  magnetic  field  B  must  be 
maintained.  Indeed,  in  Example  2  of  Chapter  5,  we  exhibit  the  effect 
on  drift  surfaces  of  not  preserving  the  divergence-free  condition  on 
solenoidal  vectors.  Also,  Anderson  et  al.  [1]  have  demonstrated  the 
rapid  loss  of  the  adiabatic  invariance  of  the  magnetic  moment  when 
V«_B    =   0  is    not    satisfied    in   their   orbit  code. 

The  divergence-free  condition  on  B  is  important  in  considerations 
about    the  stability    of    flux    surfaces,    and    this    may    be   understood    within 
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a  general  framework  in  the  theory  of  Hamiltonian  systems  that  are 
nearly  integrable.  Of  particular  importance  is  the  so-called  KAM 
theory  of  Kolmogorov,  Arnold  and  Moser  which  concerns  itself,  in  part, 
with  the  question  of  the  existence  of  invariant  curves  of  certain 
dynamical  systems.  It  is  the  solenoidal  property  of  _B  that  allows  us 
to  apply  the  KAM  theory.  We  first  describe  what  an  invariant  curve  is. 
Let  us  consider  a  transformation  T  defined  on  a  planar  cross  section  C 
of  the  torus  on  which  a  field  line  trajectory  is  initiated  at  x  and 
followed  once  around  the  torus  until  it  intersects  the  plane  again  at 
T2C^.  The  iterates  T'^x  of  Tjt  obviously  describe  the  successive 
intersections  of  the  field  line  with  the  plane  C  as  it  goes  around  the 
torus.  A  curve  T  in  C  is  said  to  be  an  invariant  curve  of  the 
transformation  T  if  the  image  T£  of  any  point  £  of  F  also  lies  on  P. 
It  is  clear  that  the  intersection  of  a  toroidal  magnetic  flux  surface 
with  the  cross  sectional  plane  C  is  a  closed  invariant  curve  of  the 
transformation  T  defined  by  the  field  line  trajectories.  Taking  as 
origin  an  arbitrary  point  0  lying  in  the  interior  of  the  area 
surrounded  by  a  closed  invariant  curve  F,  we  may  measure  the  angular 
variation  6  ^^  between  two  successive  iterates  T'^"'"  x  and  T'^x  on  T.  If 
the 


1  ""^ 

lim     1  y      9, 

n  ' 

n+oo  m=0 


exists,  then  its  value  a    is  independent  of  the  choice  of  0;  in  short,  a 

=  a(r).   a(r)  is  called  the  rotation  number  of  the  invariant  curve  T. 

The   transformation  T  defined  by  the  field  line  trajectories  has 

the  important  property  that  it  is  measure-preserving,  where  the  measure 
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is  the  magnetic  fluK.  To  be  precise,  let  y  be  any  closed  curve  on  the 
plane  C  which  bounds  an  area  Z  of  C.  Let  Ty  denote  the  closed  curve  in 
C  which  is  the  image  of  y  under  the  transformation  T;  and  let  us  write 
the  area  of  C  spanned  by  Ty  as  TE .  Our  claim  that  T  is 
flux-preserving,    therefore,    asserts    that 


JJ    B^  •    n   da      =     J  J    ^  *    n   da 
Y.  TZ 


where  n  is  a  unit  normal  to  the  plane  C.  The  proof  is  simple:  note  that 
Z  and  TZ  are  merely  the  transverse  end  faces  of  a  tubular  volume  V  in 
the  torus  whose  lateral  boundary  is  covered  by  field  lines.  The 
divergence  theorem  applied  to  the  volume  V  in  the  form 


0  =  J  J  j  V  •  B^  dx  =  -  jj    B^-  n  do  +  j  J  1*  n  da 
V  Z  TZ 


proves  the  claim  that  T  is  flux-conserving.  Since  the  magnetic  fields 
we  are  concerned  with  have  a  strong  toroidal  component,  we  may  assume 
that  _B*n_is  of  one  sign.  This  allows  us  to  define  an  "area"  element  dA 
by 

dA  =  B_  •  n^  do  . 

In  this  sense,  we  may  refer  to  T  as  an  area-preserving  transformation. 
We  wish  to  emphasize  that  T  is  area-preserving  as  a  result  of  the 
solenoidal   character   of  ^. 

The  connection  between  invariant  curves  of  T  and  its 
area-preserving  property  is  elucidated  by  a  theorem  of  Moser  [17]  and 
will    be    explained   briefly. 
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Suppose      the      area-preserving      transformation  Tq   of    the    plane    into 
itself    possesses   a  nested    set    of      closed      invariant      curves      T      In      the 
annular      region     bounded     on  the    inside   by  an  invariant    curve    F      and  on 

a. 

the    outside   by  an  invariant    curve    T,     .      Furthermore,    suppose      that      the 
rotation  numbers,   a  (F  ) ,    associated  with   these   curves    satisfy 


0  <   a(F^)    <  c.(F)    <  a(F^,) 


and   are  monotone.    If   e   is   a   small   positive   number   such   that 
ct(Fj^)  -  a(F^)  >  2e  ,  then  Moser  showed  [17]  that  for  any  number  w  in 


a(r^)  +  G   <  w  <  a(T^)    -   G 


which  is  not  closely  approxiraable  by  rational  numbers,  a  sufficiently 
smooth  area-preserving  transformation  T  of  the  plane  into  itself  which 
satisfies 


IIT^    -   TqII     <   G 


will   have   an   invariant   curve  F,    with  oo    as    its    rotation   number. 

The  effect  of  T  on  regions  where  the  invariant  curves  of  Tq  have 
rotation  numbers  that  are  near  rational  values  is  uncertain,  and  the 
only  definite  statement  that  can  be  made  is  that  the  widths  of  these 
zones   of    instability   are    small,   on  the   order  of   g    . 

The  field        line        trajectories        define        an        area-preserving 

transformation  T  of  the  cross  sectional  plane  C  of  the  torus  into 
itself.  As        we        liave        noted        previously,      this      is      due      to      the 

divergence-free  property    of    B.      The   KAM  theory,    and    Moser's      result      in 
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particular,  assures  us  that  If  the  unperturbed  magnetic  field  possesses 
a  nest  of  flux  surfaces,  a  slightly  perturbed  field,  which  is 
divergence-free,  will  still  have  many  closed  flux  surfaces. 
Interspersed  among  these  and  corresponding  to  regions  with  flux 
surfaces  belonging  to  the  unperturbed  field  that  have  rotation  numbers 
which  are  nearly  rational,  there  will  be  small  zones  of  instability 
where  the  flux  surfaces  of  the  perturbed  field  may  exhibit  complicated 
geometries. 

Now,  we  can  view  any  approximation  B  to  a  magnetic  field  B^  as 
being  a  perturbed  form  of  _B.  Our  goal,  of  course,  is  to  recover  as  far 
as  possible  the  flux  surfaces  of  B^,  if  these  exist,  using  our 
approximate  R  .  Obviously,  we  could  not  hope  to  do  so  exactly;  but 
could  we  hope  to  do  so  at  all?  The  KAM  theory  tells  us  that,  subject  to 
some  mild  smoothness  conditions  on  B  ,  if  we  ensure  that  B^  is 
divergence-free,  then  we  can  expect  that  B  will  have  many  flux 
surfaces  which,  moreover,  are  close  to  those  of  _B.  The  theory  tells  us 
nothing  about  what  to  expect  if  B  were  not  divergence-free.  In  other 
words,  by  maintaining  V«B  =  0,  we  give  ourselves  a  better  chance  of 
recovering  flux  surfaces  if  they  already  exist  for  B^.  If  B_  had  no  flux 
surfaces  to  begin  with,  then  we  have  lost  nothing  by  maintaining 
V.B^    =   0. 

The  KAM  theory  applies  to  drift  surfaces  as  well.  As  Boozer  has 
formilated  the  guiding  center  equations  [9]  and,  as  is  evident  from 
equation  (4.12),  the  guiding  center  equation  for  a  circulating  particle 
derives  from  a  solenoidal  vector  of  the  form  _B  +  r^K,  where  r^^  is  a 
small  nimiber.  The  KAM  theory  tells  us  that  if  a  nest  of  magnetic 
surfaces  exist,    then    there   are   nearby   drift    surfaces. 
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In  this  paper,  we  shall  develop  a  method  for  following  field  lines 
and  guiding  center  orbits  of  isolated  charged  particles  using  an 
equilibrium  magnetic  field  that  is  computed  by  the  code  of  Betancourt 
and  Garabedian  [3].  The  crux  of  the  problem  is  to  obtain  an  adequate 
approximation  to  the  magnetic  field  which  is  divergence-free  and  which 
permits  an  accurate  and  efficient  numerical  integration  of  the  ordinary 
differential  equations  that  describe  the  field  lines  and  orbits.  We 
shall   adopt   Boozer's    formulation   of   the    guiding   center   equation    [9]. 

The  Be tancourt-Garabedian  code  is  used  to  compute  an  equilibrium 
over  a  rectangular  grid  in  a  unit  cube  0<^s<^l,  0<^u<^l,  0£.v<^l 
into  which  the  plasma  or  vacuum  region  has  been  assumed  to  be  mapped. 
In  this  coordinate  system,  s  denotes  a  radial-like  variable,  and  u  and 
V  are  angular  variables  with  unit  periods.  The  transformation  between 
the  cv.ibe  and  physical  space  is  also  computed  as  part  of  the  solution  by 
the  equilibrium  code.  We  shall  describe  some  of  the  details  of  this 
code,  the  coordinate  system  (s,u,v)  and  tlie  transformation  in  the 
following  chapter. 

The  numerical  magnetic  field  obtained  from  the  equilibrium  code  is 
then  approximated  in  the  coordinates  s,u,v  by  a  least  squares  method  in 
a   space   of    divergence-free  vectors   which   have   the    form 

Vij;    X    Vu  +  Vx    X    Vv 

where  ij;  and  x  ^re  polynomials  in  s,  e^  and  e^-^^.  The  resulting 
representation  of  the  magnetic  field  R(s,u,v)  leads  to  very  simple 
systems  of  ordinary  differential  equations   for   tlie   field   lines   and 
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guiding        center        orbits.  The        approximation        procedure     and      the 

differential   equations   are    discussed    in  Chapters    3   and  4. 

We  have  implemented  our  method  of  approximating  B,  following  field 
lines  and  guiding  center  orbits  in  two  computer  programs,  and  we  have 
applied  them  to  study  flux  surfaces  and  drift  surfaces  for  stellarator 
field  configurations  computed  from  the  equilibrium  code.  Chapter  5 
presents  the  results  of  our  calculations,  and  Chapter  6  describes  the 
operation  of  our  codes.  We  have  also  included  the  FORTRAN  listings  of 
our  two  codes.  These  programs  have  been  written  to  be  fully  compatible 
with   the  Be tancourt-Garabedian    code. 
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2.    THE    COMPUTATIONAL   SCHEME  OF   BETANCOLTRT   AND  GARABEDIAN 

The  numerical  procedure  to  analyze  the  equilibrium  and  stability 
of  plasmas  that  was  developed  by  Betancourt  and  Garabedian  has  been 
documented  in  [3].  Various  refinements  and  applications  have  appeared 
in  subsequent  papers  [4,5,6,7],  In  this  chapter,  we  give  a  short 
description  of  their  procedure  since  the  method  we  have  adopted  for 
representing  the  equilibrium  field  output  by  their  code  depends  on  an 
understanding  of  their  coordinate  system.  Moreover,  the  two  FORTRAN 
programs  that  we  have  included  at  the  end  of  this  paper  to  implement 
our  method  are  fully  compatible  with  the  equilibrium  code.  It  is  not 
our  intention  to  give  a  complete  exegesis  of  the  Intricacies  and 
capabilities  of  the  Betancourt-Garabedian  method,  many  as  there  are, 
and  in  particular,  we  shall  not  discuss  the  stability  analysis  by  their 
code  of  a  computed  equilibrium,  this  being  beyond  our  purposes. 
Rather,  vre  shall  be.  content  with  presenting  so  brief  a  synopsis  of  the 
equilibrium  computation  as  will  facilitate  a  full  understanding  of  our 
ideas   and  their    subsequent   applications. 

The  method  of  the  equilibrium  code  is  based  on  a  variational 
principle  of  ideal  magnetohydrodynamics  [15]  which  characterizes 
equilibria  as   stationary  points    of    the    potential  energy 

1   .2  ^     P 


(2.1)  E  =  JjJ     [    i-B^   +-^]    dV 


subject  to  flux  constraints  to  be  described  shortly.  Here  B  is  the 
magnitude  of  the  magnetic  field  B^,  p  is  the  fluid  pressure  defined  in 
terms  of  the  density  p  by  the  equation  of  state,  p  =  p^ ,  and  y  is  the 
ratio   of    specific  heats. 
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In    the    plasma,    B^  Is   represented   as 

(2.2)  2     "     Vs  X    Vc 

where  s  and  ^  are  two  functions  on  which  flux  constraints  are  Imposed. 
The  function  s  is  assumed  to  be  single-valued  and  to  define  a  nested 
family  of  toroidal  surfaces,  s  =  constant,  for  values  of  s  ranging  from 
zero  to  one.  The  innermost  surface,  s  =  0,  degenerates  to  a  closed 
field  line,  called  the  magnetic  axis,  and  the  outermost  surface,  s  =  1, 
is  the  interface  between  the  plasma  and  the  vacuum.  Any  point  on  a 
surface  s  =  constant  is  characterized  by  two  angular  coordinates  u  and 
V  with  unit  periods,  corresponding  to  the  poloidal  (or,  short)  and  the 
toroidal  (or,  long)  directions  around  the  torus,  respectively. 
The  flux    function  ^    has    periods    given  by 

(2.3)  C    =  -  u  +  u(s)v  +  X(s,u,v) 

where  \    is  single-valued,  and  \i  ,    known  as  the  rotational  transform,  is 
defined  to  be  the  ratio 


fp(s) 
(2.4)  u(s)  =  -.U- 

t^(s) 


fp(s)  and  f'j'(s)  are  the  prescribed  poloidal  and  toroidal  nvagnetic 
fluxes  subject  to  which  the  variation  of  the  energy  E  is  made.  The 
rotational  transform  characterizes  the  twist  of  a  magnetic  field  line 
on  the  toroidal  surface  s  =  constant  and  equals  a/2Ti  ,  where  a  is  the 
rotation   number  which  we    discussed   in   connection  with    the     KAM      theory. 
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In   the   usual  notation  of  the  literature  on  plasma  physics,  q  =  \i~     is 
called  the  safety  factor.   Also,  the  rotation  number  a  is  often  written 
as   I   (iota),  but  we  have  not  used  this  symbol  to  avoid  confusion  with 
the  complex  number  i  =  /-I  which  we  shall  have  occasion  to  use. 

In  the  vacuum  region,  B^  is   represented   as   a  harmonic   vector, 
derived  from  a  scalar  potential  function  (^ , 

_B   =  V(J)  . 

The      function  (})    may   have    singularities   associated  with   current-carrying 
coils    that   generate  the    magnetic    field    [7]. 

The  constraints  on  the  variation  of  the  potential  energy  E  are 
five    [3]: 

1.  V«_B    =   0  everywhere  ■•  -  ,  .  • 

2.  The      toroidal      flux,       f^(s),    and   the    poloidal    flux,    fp(s),    of    the 
magnetic    field  B  within  each   surface    s=  constant   are    fixed. 

3.  The  mass,  M(s),    within  each   surface    s   =  constant    is   also    fixed. 

4.  The    total    toroidal   and    poloidal    fluxes    in   the    vacuum     region      are 
fixed . 

5.  The      free    surface    between   the    plasma   and   vacuum  regions    is    a   flux 
surface .  . 

Under  these  conditions  it  can  be  shown  that  the  Euler  equations  of 
the  variational  problem  reduce  to  the  magnetos tatic  equations,  and  that 
the  energy  E  is  minimized  if  p  and  p  are  functions  of  s  alone  [3].  In 
the   vacuim    region,   we  will   have 

A<t.    =    0  '■ 
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and  in    the    plasma    region, 

(2.5)  ■  Vp   =    J  X    B 

1      o 
where     J_  =  V   x    B_.      Moreover,    the   quantity  _  B      +  p    is   continuous   across 

the    free    boundary   between    the    plasma   and   vacuum  regions.      Observe      that 

equation   (2.2)    implies    that    the    last    magnetos  tat ic   equation,    namely, 

(2.6)  V    •    B_     =     0 

is  automatically  satisfied.  The  fact  that  the  pressure  and  density 
functions  depend  on  s  only  means  that,  without  loss  of  generality,  we 
may  include  this  requirement  as  an  additional  constraint  on  the 
variational  problem.  This  has  an  important  consequence  which  we  will 
discuss    in  due  course. 

If      ('^i)X2,Xo)   are  rectangular  Cartesian    coordinates    in    the    plasma 
and 


3 (x^ ,X2,X3) 
8 (s ,    u,   v) 


8  (s,   C  ,    X.) 
D.    =        ,  \ 

J         8 (s ,    u,    v) 


then  equation  (2.1)  for  the  energy  can  be  reformulated  as 


(2.7)   E  =  jjj  ^^      ^   ds  du  dv  +  -L  J  __ii_Ll2 

2n  ^^  N  JJ  D  du  dv)Y-l 


The  crux  of  the  numerical  method  is  a  finite  element  quadrature 
approximation  to  this  last  equation  for  the  energy  based  on  a 
rectangular  mesh  over  the  cube  0<s<I,  0<u<l,0<  v  <  1.  The 
presence     of       the      .Tacobians      T)      and   n      i  ,i   thi^    equation    means    that    the 
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t rans formation   between   the    flux      coordinates      s,u,v     and      the      physical 
coordinates    x, jXojXo    has  to  be   computed   as  part    of    the    solution. 

Let  (Rq.Q.z)  be  the  cylindrical  coordinates  of  the  point 
(x^,X9,X3)  in  the  plasma.  If  A  denotes  the  large  radius  of  the  torus, 
then 

Xj  =  (A+  r)  cos  e 
X2  =  ( A  +  r)  s  in  Q 
X-j    =    z 

where  r  =  Rq  -  A.  Ostensibly,  the  energy  E  is  then  a  functional  of  the 
four  variables  ^,r,z,9  which  agrees  with  the  fact  that  the  Euler 
equations  (2.5)  and  (2.6)  for  the  plasma  region  comprise  a  system  of 
four  partial  differential  equations  for  the  function  p  and  the  three 
components  of  B.  The  choice  of  B  as  given  by  equation  (2.2)  and  the 
additional  constraint,  p  =  p(s),  serve  to  reduce  the  order  of  this 
system   by   two   since    B  is  automatically  divergence-free  and 

(2.8)  B^  •    Vp     =      0 

holds.  More  significantly,  relation  (2.8)  implies  that  the  real 
characteristics  of  the  Euler  equations  have  been  eliminated,  and  what 
remains  is  a  system  of  two  elliptic  partial  differential  equations  for 
which      well-posed      boundary      value      problems      may      be        set.  Another 

consequence  of  the  reduction  of  the  order  of  the  Euler  equations  is 
that    now  the    function  9    may   be   chosen   to   satisfy 
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(2.9)  e    =   2irv 

thus   diminishing   the   number   of   unknown    functions    by   one. 

There  are  three  reasons  why  it  would  not  be  advantageous  to  treat 
the  variable  c,  as  a  given  function  of  s,u  and  v,  and  to  minimize  the 
potential  energy  over  all  periodic  mapping  functions  r  and  z.  First  of 
all,  the  flux  function  c,  is  not  unique  since,  as  equation  (2.2)  shows, 
an  arbitrary  function  of  s  may  be  added  to  it  without  affecting  B^. 
Second,  the  magnetic  axis  is  a  singular  curve  of  the  coordinate  system, 
and  it  would  be  difficult  to  write  accurate  difference  equations  for  r 
and  z  in  its  vicinity.  Last  of  all,  the  boundary  conditions  for  r  and 
z   are   nonlinear. 

A  simpler  formulation  is  achieved  by  minimizing  the  potential 
energy  in  the  plasma  region  as  a  functional  of  the  variables  ^  and  a 
dimensionless  radial-like  variable  R,  defined  by  means  of  the  following 
f  omiilae: 

(2.10)  r      =   rQ(v)   +  R(s  ,u,v)  [  r^  (u  ,v)    -  ro(v)  ] 

(2.11)  z      =  zo(v)   +  R(s,u,v)[zi(u,v)    -  zo(v)] 

Here,  r  =  rg(v)  and  z  =  zq(v)  are  the  equations  of  the  magnetic  axis, 
and  r  =  r^(u,v)  and  z  =  zj(u,v)  are  the  equations  of  the  free  boundary. 
The  function  R(s,u,v),  which  is  periodic  in  u  and  v,  equals  one  when 
s  =  1,  and  is  zero  when  s  =  0.  For  cases  with  no  vacuum  region,  the 
functions  r  =  r^(u,v)  and  z  =  zj(u,v)  become  equations  for  the  outer 
conducting  wall  and  must  be  input  to  the  equilibrium  code.  The  energy 
E,      as      given  by  equation    (2.7)is    then    considered    to  be   a    functional   of 
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the  variables   c,  ,^,Tq^zq,    and    its  rainimura  is  sought   by   defining   paths    of 
steepest    descent    and   employing    the    conjugate    gradient    method. 

^«Jhen  a  vacuum  region  is  present,  the  location  of  the  free  boundary 
has  to  be  found  through  the  contribution  of  the  vacuum  region  to  the 
value   of    the   energy,    and    the    free    boundary   condition 


1      ? 
—  B     +  p   =  constant. 


We  shall  not  describe  how  this  is  carried  out  since  we  are  primarily 
interested  in   following    field   lines    and   orbits    in    the    plasma    region. 

At    the    termination  of  a   successful  run   of    the   equilibrium  code    for 

a      case     when    the    vacuum  region   is    absent,    an  equilibrium  is   assumed    to 

have   been  achieved,   and    the    values    of    the    minimizing    functions   c, ,    R,    r^ 

and  Zq  are,  therefore,  known  at  equally  spaced  grid  points  in  the  unit 
(s,u,v)  cube.  These  data  constitute  the  raw  materials  on  which  our 
program  will  act. 

It  is  appropriate  to  emphasize  some  of  the  features  of  the 
Re  tancourt-Garabedian  method.  It  tias  recast  the  variational  principle 
in  such  a  way  as  will  lead  to  nontrivial  well-posed  problems  in 
toroidal  geometry  with  no  spatial  symmetries.  This  fact  and  the  simple 
way  in  which  the  constraints  and  boundary  conditions  have  been 
incorporated  into  the  variational  principle  are  very  important  for 
sound  numerical  analysis  and  computing.  The  generally  good  agreement 
between  results  computed  by  the  equilibrium  code  and  tVieoretical 
predictions  as  well  as  experimental  data  suggests  that  such  concerns 
may  not  be  entirely  without  pne  ri  t .  The  significance  of  the  assumption 
of      nested       flux    siirfaces    in    their    formulation    vis-a-vis    the  KAM   theory 
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has  been  discussed  at  greater  length  by  Garabedian  in  [10]  and  also  in 
[7].  The  philosophy  may  be  simply  stated:  even  though  there  may  be 
regions  of  ergodicity  of  the  field  lines  where  no  flux  surfaces  are 
swept  out,  the  KAM  theory  predicts  that,  under  favorable  conditions,  an 
infinite  number  of  nested  toroidal  flux  surfaces  will  still  remain,  and 
will  predominate  over  the  smaller  regions  of  ergodic  behavior.  The 
finite  nunber  of  nested  flux  surfaces  implicit  in  the  computational 
grid  used  by  tlie  equilibrium  code  may  be  thought  of  as  constituting  a 
selection  from  this  infinite  set  of  nested  surfaces.  Nevertheless,  it 
is  still  necessary  to  have  independent  ways  of  testing  for  the  break-up 
of  surfaces;  these  have  been  developed  and  incorporated  into  the 
equilibrium  code  [7].  In  particular,  the  method  of  following  field 
lines  to  be  described  in  the  next  chapter  leads  to  an  additional 
quantitative   assessment    of    the   adequacy  of    the   equilibrium  code. 
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3.  A  METHOD  TO  STUDY  MAGNETIC  SURFACES  IN  TOROIDAL  PLASMAS 

BY  FOLLOWING  FIELD  LINES, 
(a)  The  Representation  of  the  Magnetic  Field. 

The  choice  of  a  representation  of  the  magnetic  field  is  governed 
largely  by  its  appropriateness  to  the  problem  at  hand.  For  example,  in 
studying  the  local  structure  of  the  magnetic  field,  it  is  often  useful 
to  write  _B  as 

B  =  VF  X  VG 

where  F  and  G  are  functions  of  position.  This  representation  of  B, 
sometimes  called  its  Clebsch  representation,  is  always  possible  in  a 
localized  region  of  the  torus.  The  standard  argument  begins  by 
defining  two  arbitrary  functions  F  and  G  on  a  planar  cross  section  C  of 
the  torus  so  that  their  level  curves  form  a  valid  coordinate  system  on 
C.  The  normal  component  of  B  is  supposed  never  to  vanish  on  this  plane 
or,  at  least,  never  over  a  localized  area  on  C.  We  then  define  the 
functions  F  and  G  at  points  on  any  other  planar  cross  section  C'  nearby 
to  C  by  using  the  field  lines  to  convect  the  values  of  F  and  G  onto 
points  of  C.  Ihe  level  surfaces  of  the  extended  functions  F  and  G 
intersect  along  field  lines  as  a  consquence,  and  this  entails  that  the 
vectors  _B  and  VF  x    VG  be    parallel.       So,    we   may   write 

_B    =  f    VF  X    VG 

for  some  function  f  of  position.  In  fact,  f  is  constant  along  a  curve 
of  intersection  of  two  level  surfaces  of  F  and  G,  that  is,  along  a 
field  line.   Indeed,  the  solenoidal  property  of  B  leads  to 
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Vf    •    VF  X    VG   =    0 

which  means    Chat    f  =   f(F,G).      We    can   then   define  F,    by 

F 
F^(F,G)       =     J     f(e,G)    dC 

to  arrive   at 


B      =     VF^   X    VG 


This  argument  applies  In  a  localized  region  of  the  torus,  and  may 
suffice  for  some  purposes.  Whether  it  is  possible  to  obtain  a  global 
Clebsch  representation  of  B  by  this  means  is  a  difficult  question  to 
answer  [2].  If  a  nested  family  of  flux  surfaces,  given  by  s  = 
constant,  exists,  then  we  may  use  these  to  define  one  of  the  coordinate 
functions,  say  F.  We  are  then  led  to  the  representation  (2.2)  of  _B  as 
used  by  the  Be tancoiirt-Garabedian  method.  Our  purpose,  being  to  find  a 
global  approximation  to  a  magnetic  field  with  no  prior  assumption  about 
the  geometry  of  the  flux  surfaces,  necessitates  a  different  form  of  a 
Clebsch  representation  for  3  output  by  the  equilibrium  code. 
Hopefully,  where  the  surfaces  are  indeed  nested,  ours  will  be  a 
faithful  representation  of  the  computed  magnetic  field,  and  where  the 
hypothesis  of  nested  flux  surfaces  fails,  our  method  may  reveal  the 
extent    to  which   this    is   untenable. 

Our  procedure  for  representing  the  magnetic  field  will  now  be 
described.  It  is  assumed  that  a  numerical  solution  of  the  variational 
problem  is  known  from  the  equilibrium  code  on  a  rectangular  grid  over 
the   unit  cube   whose   coordiates   are    s,    u   and      v,      as      described      in      the 
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previous  chapter.  In  particular,  the  values  of  the  flux  function  c, 
appearing  in  equation  (2.2)  have  been  conputed  at  the  grid  points.  We 
may  then  use  these  values  to  compute  finite  difference  approximations 
to  the  two  components  D_B«Vu  and  DB^-Vv  of  the  magnetic  field  B^  at  grid 
points    by  means    of    the    fomulae 


(3.1)  DB   •    Vu  =    -^ 

-  3v 


and 


(3.2)  DB  .    Vv  =       ii 

8u 


v-7hlch   follow  from   equation   (2.2).      4ere  the   Jacobian   D   is 


8(x,,X9,X3) 

(3.3)  0  =        ^  _     =      (Vs   •    Vu  X   Vv)    ^ 

3 (s ,    u,   v) 


We  recall  that  whereas  the  function  r,  ,  as  expressed  by  equation  (2.3), 
is  mi.ilt  ivalued  in  the  poloidal  u  and  toroidal  v  angular  variables,  the 
magnetic    field   is  periodic    in   these    variables    with   unit    periods. 

Our  method  is  based  on  the  idea  of  approximating  the  magnetic 
field  computed  by  the  equilibrium  code  in  a  space  of  divergence-free 
vectors.  The  technique  may  be  considered  as  an  application  of  the 
method  of  weighted  residuals  by  point  collocation,  but  it  will  be 
apparent  that,  in  our  case,  this  reduces  simply  to  the  least-squares 
fitting   of    data    by   polynomials  and    trigonometric    functions. 

The  central    idea   in   our    representation    is    the   expression   of    B  as 
B_   =  V\J;    X    Vu  +  Vx    X    Vv 
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A  simple  argument  will  show  that  we  can  always  do   this.    Indeed,   the 
solenoidal  magnetic  field  B_  is  derivable  from  a  vector  potential  A_  [16] 
as 

B  =  V  X  A 
We  may  express  the  vector  A_  in  covariant  form  by 

A  =  A^Vs  +  A2VU  +  A3VV 
since  Vs,  Vu  and  Vv  are  linearly   independent.    If   we   now   define   a 
scalar  function  (J)  by  means  of 

s 

(j)(s,u,v)   =  j   A^(t,u,v)  dt 


then    it    follows    immediately    that 

AjVs      =      V(ti    -  4)^Vu   -   (Jj^Vv 
where     <p ^  =  _I-  and  <^ ^   =  —I-   .      Consequently,    by    choosing  41   =  A9    -  (^^   and 
X   =  A^   -  ({)„    ,    we   will   arrive   at    the   desired   expression 


(3.4)  B^     =     Vij)   X    Vu  +  Vx   X   Vv 


We  shall  assume  that  ^    and  x  have  the  following  expansions 


(3.5)     :i,    =1   a^^^Pj^(s)  ei2^(rau-"^) 
(3-^)     X  =1    h:nn%^^^    ei2.(im-nv) 


where  the  functions  P^(s)  and  Qj^  (s )  are  polynomials,  and  the  olq^^  and 
^£mn  ^^^  constants.  The  equations  (3.4),  (3.5),  (3.6)  describe  our 
representation  of  the  magnetic  field.  To  motivate  our  choices  for  \|; 
and  x>  ^6  observe  that  the  forraLilae  for  the  two  components  DB'Vu  and 
DB_«Vv  of  _B  which  follow  from  the  last  three  equations  are 
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(3.7)  DB-Vu  =  -I   ej^ninQ^(s)  e^^irCmu-nv) 
and 

(3.8)  DB.Vv  =   I   a^^^P^(s)  e^^Tr  (tmi-nv) 

where  the  primes  ( ' )  denote  differentiation  with  respect  to  the 
indicated  variable.  Now,  the  fact  that  a  sufficiently  smooth  magnetic 
field  B^(s,u,v),  which  is  periodic  in  the  variables  u  and  v,  will  have 
regular  components  DB«Vu  and  D_B«Vv  that  allow  Fourier  expansions  as 
those  of  equations  (3.7)  and  (3.8)  leads  us  to  believe  that  we  have 
chosen  suitable  formulae  for  the  functions  \\)    and  x» 

The  procedure  for  approximating  the  magnetic  field  computed  by  the 

equilibrium  code  consists  of  determining  the  constants   a„     and  6„ 

^  °  xmn  "^xmn 

from     the      data   points.      Tt    follows    from  equations   (3.2)    and    (3.8)    that 

the   coefficients     ct^^^^^     may      be      found      by      Fourier     analyzing      in      the 

variables      u     and      v,      and      polynomial      least-squares      fitting   in  s   the 

equally  spaced   data    points    of    the    function  Jz.  .      The    coefficients      3o„„ 

8u  x.ran 

are  determined  by  the  same  technique  applied  to  the  data  points  of  the 
function  _ —  .  These  constants  a^^^  and  S>ji^^  depend  on  the  particular 
grid  size  at  which  the  equilibrium  was  calculated.  The  successive 
determination  of  these  constants  from  equilibrium  fields  computed  over 
various     meshes      allows     us   to  extrapolate   their   values   to   the    limit   of 

zero  mesh  size.      Once  we  have  computed    the    coefficients   a„        and  6n 

Jimn     '"^K.mn  ' 

either  with  or  without  extrapolation,  the  expressions  (3.7)  and  (3.8) 
for  DB^'Vu  and  DB»Vv  are  known.  The  simplest  way  then  to  obtain  a 
formula  for  DB^»Vs  is  to  integrate  with  respect  to  the  variable  s  the 
equation  V«B  =  0  written  as 
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(3.9)  i_DB«Vs      =      -  _L  DB^Vu   -  —  DB-Vv 

dS       —  dU—  dV 


Carrying   out    this    straightforward   calculation,    we   arrive  at 

(3.10)      DB-Vs  =   f(u.v)    +1   ^s,^n^?iis)-?^isQ)]e^^''^'^-''''^ 

^l    WQ.(^)-Q.(so)]e^2.(™-nv) 

where 

"Jlmn     =      i2Tin  a^^^^ 

^ilmn     =      i2Trm  B^^^ 

and  f(u,v)  is  any  periodic  function  of  the  variables  u  and  v.  If  we 
assume  that  the  surface  s  =  Sq  ,  for  some  particular  value  of  Sq  ,  is  a 
flux  surface,  so  that  DB^«Vs  =  0  on  that  surface,  then  f  is  identically 
zero,    and,    therefore,    the  component   DB^«Vs    is   uniquely   determined. 

The  fact  that  we  have  been  able  to  obtain  the  formulae  for  the  two 
components  DB^»Vu  and  DB«Vv  independently  of  each  other  is  a  direct 
consequence  of  the  way  we  have  chosen  to  represent  B^  through  equations 
(3.4),  (3.5)  and  (3.6).  The  divergence-free  condition  on  B^  serves  to 
determine  the  remaining  component  D_B«Vs  up  to  an  arbitrary  periodic 
function  in  the  variables  u  and  v;  this  is  also  reflected  in  the  fact 
that  such  a  function  may  be  added  to  either  i|;  or  x  without  changing  the 
formulae  (3.7)  and  (3.8)  for  the  component  D_B«Vu  and  D_B»Vv.  To 
determine  DB«Vs    uniquely,    we   assume    the   existence   of    one   flux   surface. 
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(b)    The  Field  Line  Equations 

If  we  assume  that  we  now  have  a  representation  of  the  magnetic 
field  described  by  (3.4),  (3.5)  and  (3.6),  then  it  follows  that  the 
system   of  ordinary  diffrential  equations 


dx 

T=     =     B(x) 


dt 


for  the    field   lines    is    given  in   the    coordinates   s,u,v  by 


ds 

-^     =      DB   .    Vs 
dt  - 


i^     =      DB   .    Vu 
dt  - 


dv 

_     =     DB  •    Vv 

dt  - 


As  is  evident  from  equations  (3.7),  (3.8)  and  (3.10),  the  right  members 
of  this  last  system  of  differential  equations  consist  essentially  of 
Fourier  series,  and  their  evaluation  is  quite  simple.  Moreover,  the 
integration  of  a  field  line  in  these  coordinates  is  very  efficient, 
since  experience  has  shown  that  we  need  retain  relatively  few  harmonics 
in  the  series.  In  practice,  we  have  successfully  employed  a 
fourth-order  Adams-Moulton  predictor-corrector  multistep  integration 
scheme  coupled  to  a  fourth-order  Runge-Kutta  starter.  To  integrate  a 
field  line  once  around  a  torus  with  inverse  aspect  ratio  (average  minor 
plasma  radius  to  average  major  plasma  radius)  of  0.1  requires  on  an 
average   10   seconds   on   the   CDC    6600  computer   at   New  York   University. 

For     purposes      such      as      plotting      the      solution      of   the   numerical 
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integration  in  physical  space,  we  must  use  the  transformation  described 
in  Chapter  2  by  equations  (2.9),  (2.10)  and  (2.11).  The  functions 
rr.{v) ,  Zn(v)  and  R(s,u,v)  which  appear  in  those  equations  are  computed 
by  the  equilibrium  co(ie  on  a  fixed  grid,  and  it  is  necessary  to  be  able 
to  approximate  their  values  at  points  not  on  such  a  grid.  Since  these 
functions  are  periodic  in  the  variables  u  and  v  with  period  unity,  we 
expand    thera   in  Fourier    series   as 

(3.11)  ro(v)      =      I   a,ei2--v 

(3.12)  Zo(v)      =      I   b^e^2^"v 

(3.13)  R(s,u,v)      =     I    c.^^s^    ei2.(mu-nv) 

and  determine  the  coefficients  a^  ,  b^^  ,  (^^^^  from  the  data  by 
least-squares    fitting   as   before. 

The  method  of  following  field  lines  described  here  has  been 
applied  to  field  configurations  in  stellarator  geometry  computed  by  the 
equilibrium  code,    and    we   refer    to   Example    1    of  Chapter    5. 

One  issue  that  we  must  address  concerns  the  question  of  whether  or 
not  the  assumption  of  nested  flux  surfaces  in  the  method  of  the 
equilibrium  code  predisposes  our  procedure  merely  to  reproduce  those 
nested  surfaces.  We  immediately  concede  that  the  direct  application  of 
our  method  to  a  magnetic  field  computed  on  a  fixed  grid  by  the 
equilibrium  code  has  never  produced  flux  surfaces  with  geometries 
different  from  those  of  this  code.  However,  this  was  not  the  case  when 
we  followed  field  lines  using  the  magnetic  field  extrapolated  from 
those  produced  by  several  runs  of  the  equilibrium  code  made  on  various 
grids.      In    that    instance,    as   is    discussed    in  Example    1    of  Chapter    5,    we 
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found      a      chain      of     magnetic      islands,    a  structure  which   is    apparently 
precluded     by      the      equilibrium        code.  This        suggests         that        the 

extrapolation  of  results  from  the  equilibrium  and  stability  code  may  be 
markedly  different  from  those  obtained  at  a  fixed  mesh.  Our  example 
exhibits  the  fact  that  we  have  not  assumed  the  nestedness  of  the  flux 
surfaces  in  our  representation.  Although  this  is  very  important,  in 
and  of  itself,  it  would  not  explain  why  our  scheme  allows  us  to  capture 
these  island  structures,  starting  as  it  does  with  a  magnetic  field  from 
the  equilibrium  code  which  does  assume  the  hypothesis  of  nested  flux 
surfaces.  We  believe  that  our  ability  to  reproduce  island  structures, 
if  these  exist,  may  be  due  to  the  smoothing  of  the  equilibrium  magnetic 
field  in  its  radial  dependence  on  s  by  our  polynomials  ?.  and  Qp  .  The 
results  of  recent  studies  by  Betancourt  and  Garabedian  [7]  have  led 
them  to  believe  that  in  their  variational  formulation  the  hypothesis  of 
nested  flux  surfaces  and  the  regularity  of  the  rotational  transform  \i 
both  cause  the  current  density  J_  =  VxB^to  be  discontinuous  in  the 
vicinity  of  a  rational  magnetic  surface,  on  which  the  rotation  number 
is  a  rational  multiple  of  2tt  .  We  recall  from  our  discussion  of  the  KAM 
theory  that  this  is  a  region  where  surfaces  may  be  most  vulnerable  to 
small  perturbations.  In  our  method,  by  approximating  the  radial 
dependence  of  all  functions  by  polynomials,  we  are  assuming  that  B^  is 
dif f erentiable  in  s  everywhere  and,  moreover,  that  J_  is  also  regular. 
This  coupled  with  the  relaxation  of  the  nested  surfaces  constraint 
allows  the  rotational  transform  to  be  discontinuous  and  the  islands  to 
appear. 
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4.    THE    A.PPLICATION   OF  THE    METHOD   TO   STUDY    DRIFT   SURFACES 
(a)   The  Guiding  Center  Equations 

The     exact      motion      of      a      charged    particle   in  an  external,    static 
magnetic   field  B   is  determined  by  Newton's    second   law   of   motion, 


dv 


m  —     =      qv   X    B 
dt  -       - 


Vtere,  v^  is  the  velocity  of  the  particle,  whose  mass  and  charge  are 
denoted  by  m  and  q,  respectively.  The  right  member  of  this  equation  is 
simply  the  Lorentz  force  on  a  charged  particle  due  to  the  magnetic 
field.  ^/hen  the  field  is  constant,  this  equation  can  be  integrated 
explicitly.  The  motion  is  found  to  be  a  superposition  of  a  circular 
gyration  about  a  field  line  with  gyrof requency  9,  =  qB/m,  and  a  uniform 
translation  in  the  direction  of  _B.  The  orbit  is,  therefore,  helical. 
The  radius  of  the  circle  of  gyration  is  called  the  gyroradius,  and  the 
center   of   gyration   is    the    guiding   center. 

When  the  field  is  inhomogeneous  but  changes  by  small  amounts  over 
distances  comparable  to  the  gyroradius,  we  can  study  the  motion  of  the 
particle  by  means  of  the  guiding  center  approximation,  detailed 
derivations  of  which  will  he  found  in  the  literature  [16,18].  This 
approximation  is  a  small  parameter  expansion  of  the  particle  orbit 
about  a  field  line.  The  small  parameter  is  the  ratio  of  the  gyroradius 
to  a  characteristic  length  (usually  the  minor  radius  of  the  plasma) 
over  which  the  magnetic  field  changes  appreciably.  The  guiding  center 
theory  describes  the  motion  of  the  instantaneous  center  of  gyration  of 
the  particle;  it  consists  of  a  component  parallel  to  the  magnetic 
field,    and  a   drift    component    perpendicular   to    the    field.      Whereas,    in  a 
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uniform  field,  the  guiding  center  orbit  coincides  with  the  field  line 
so  that  the  drift  surfaces  are  identical  to  the  flux  surfaces,  in  an 
inhoraogeneous  field,  the  guiding  center  will  drift  across  field  lines 
causing  the  drift  surfaces,  if  they  exist,  to  deviate  from  the  flux 
surfaces.  It  is  this  slower  drift  of  the  orbit  that  may  prevent  the 
long   term  confinement    of   isolated  charged    particles. 

In  this  work,  we  have  investigated  the  nature  of  drift  surfaces  by 
integrating  the  first  order  system  of  ordinary  differential  equations 
which  describes  the  trajectory  of  the  guiding  center.  We  have  followed 
Boozer's  formulation  of  this  system  [9],  since  it  lends  itself  readily 
to  our  technique  of  representing  divergence-free  vectors  in  Clebsch 
form.       Tfeozer's    equation    for   the    guiding   center    is 


dx  ^" 


(4.1)  ^     =     ^[B   -fV   X    p,|Bj 


where 


X  =  the  position  vector  of  the  guiding  center 

B  =  the  magnetic  field  vector 

B  =  111 

V||  =  the  magnitude   of    the    parallel    (to  _B)  component    of 

the    particle's   velocity 
mvii 


^"  qB 

m        =     the  mass   of   the    particle 

q        =      the  charge  of   the    particle 


Boozer's      equation      (4.1)      differs      from   the   usual   equation  of   the 
guiding   center    [16]    in  its    expression    for   that    component    of    the    guiding 
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center      velocity     v^  which   is   parallel    to  _B.      Both  formulations   agree    on 
the   expression    for   the    drift    velocity   v,     ,    which   is    the   component    of      v 
perpendicular   to  _B.      The  usual   expression    for   the    parallel  component    of 

V|| 

_v  is  B   .      Observe   that   this    is    a   quantity   of    the    zeroth   order    in   the 

B 

gyroradius.    By   adding   an  appropriate  term  of  the  first  order  in  the 

V|| 

gyroradius   to  —  B   ,    Boozer  was  able    to  arrive   at   his    equation   (4.1). 
B 

1  7  1  "^ 

The  energy  W   =  _  mvf    +—  mvT  of    the   particle    is   conserved,   and    the 

1         7 
magnetic   moment   P    =  -=-  mvr/B    is   an  adiabatic   invariant    of    the   motion  and 

so  also  constant.      The  two   constants    of    the   motion,   W  and   y,      serve      to 

determine  v..    as   a    function  of    position,   since 


(4.2)  V||    =  ±    [|  (W  -  yB)]^/2 


Observe  that  we  have  used  the  same  symbol  u  to  denote  both  the 
magnetic  moment  and  the  rotational  transform.  No  confusion  will  arise 
as  the  meaning  of  each  occurrence  of  y  will  be  apparent  from  the 
context,    or   explicitly  noted. 

For  given  values  of  the  energy  and  the  magnetic  moment  there  are 
roughly    two    categories    of    orbits: 

1.  In  regions  of  high  field  strength,  v..  may  vanish  and  the  particle 
will  be  reflected  back  along  a  field  line.  This  behavior  is  the 
principle  behind  particle  confinement  in  mirror  systems,  but  noay 
also  be  present  in  toroidal  systems  due  to  the  variation  of  B. 
The  particle  is  said  to  be  in  a  trapped  state  when  it  bounces 
back  and    forth   between  reflection   points    where   W   =  yB. 

2,  The   magnetic  moment  may  be  too   small   to  allow      the      vanishing      of 
V|.  ,      that    is,    W  >  yB  along   the   orbit.      In    this   case,    the    particle 
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will  traverse    r.he    torus  unimpeded  by  reflections,    and    it   is      then 
said   to  be    untrapped,   or  is    called   a   circulating   particle. 
These      two      categories      are   not    nutually  exclusive,    and   a  particle 
may   spend  part    of    its    time    in   one    or    the      other      state.        For     a     given 
field      geometry,      the      guiding      center      orbit     may   be  quite    complicated 
depending  on   the    two  parameters   W  and   y. 

The  possibiity  that  v..  may  vanish  at  points  along  the  guiding 
center  orbit  would  seem  to  cause  difficulty  in  the  integration  of 
equation  (4.1),  especially  because  of  the  singular  nature  of  the 
derivative  of  p..  at  such  points.  This  singularity  is  nonessential, 
however,   for  we   can  rewrite    (4.1)   as 


dx 


2 
At      reflection      points      of      the      orbit,      the      derivative        of        P||         is 

well-behaved. 

In     addition,      there      is    a   practical    question  of  how  to  choose   the 

proper   sign  of   the    square    root   in    the   expression   (4.2)      for     v.,       during 

the    integration.      As  Boozer  has    suggested    [9],    this    can  be   conveniently 

resolved      by      integrating,      in     addition      to      tlie        system        (4.1),        a 

differential   equation   for  p..     ,    namely, 


(^•^)    ,  ^  =  ^[1-    V    ^P?      +     P„    V    X    B..    Vl-pif] 


1      2 
An      expression   for   the   vector  V  —  p ,,    ,    appearing  on   the    right-hand  side 

of    (4.4),   may  be   derived    from    (4.2)as 
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_    1      2  m     r    uB-2W    T 


It  should  be  noted  that  the  numerical  integration  of  the  augmented 
system  of  differential  equations  (4.3)  and  (4,4)  can  no  longer  be 
expected  to  conserve  the  energy  W.  Fortunately,  we  can  use  this  fact  to 
advantage,  since  the  energy  may  be  calculated  at  each  step  of  the 
integration  of  (4,3)  and  (4.4)  from  the  solutions  2i(t )  and  P||(t)  ,  and 
its  deviation  from  constancy  will  provide  a  useful  measure  by  which  to 
assess  the  accuracy  of  the  numerical  scheme.  We  are  grateful  to 
Dr.    Boozer   and   Dr.    Kuo-Petravic    for   their   advice   on    this    point. 

When  a  particle  is  circulating,  so  that  v.,  does  not  vanish  at  any 
point  along  the  orbit  of  its  guiding  center,  the  factor  v., /B  in 
equation  (4.1)  is  of  no  consequence,  and  we  may  regard  the  guiding 
center  orbit  of  such  a  particle  as  being  a  field  line  trajectory 
belonging  to  the  solenoidal  vector  field  B_  +  V  x  p  g  .  This  vector  is 
a  modification  of  _3  by  the  terra  V  x  P||B^  ,  which  as  we  shall  see  in 
section  4b,  is  of  the  order  of  the  gyroradius.  The  KAM  theory 
suggests,  therefore,  tliat  if  magnetic  surfaces  exist,  and  if  the 
gyroradius  is  sufficiently  small,  then  there  will  be  nearby  drift 
surfaces  belonging    to   the    divergence-free  vector  B^  +  V   x   p    b   . 

The  significance  of  Boozer's  formulation  to  our  work  lies  in  the 
fact  that  the  right  member  of  the  guiding  center  equation  ((4.1)  or 
(.4.3)  )is  expressed  as  a  sum  of  terms  that  involve  solenoidal  vectors. 
This  allows  us  to  write  each  of  these  vectors  in  our  Clebsch 
representation. 
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1      ? 
For  example,    let  us    denote   the    vector  ^   ^   -y  Pif  B^  by   H,      and     by     H 

its  Clebsch    represent.^tlon,   that    Is 


where 


H    =  Vi|;  j^  X   Vu     +     Vx^   X   Vv 


,,,    _  Y    r  P    (<i)    „i2TT(mu-nv) 

1^1-  i    *^mn  ^i'^^^    ^ 


and 


V  r\    f    \      i2TT  (mu-nv) 

Just  as  in  the  case  of  _B,  the  coefficients  Y^^^^  ^'^'^  ^SLmn  ^^  '-^^ 
fomulae 

'^ii-Vu=    -I   Y,,„Q^(s)    e^2-(  — ) 
Dii-Vv=       I    5,,„„P^(s)    ei2-(  — > 

are  determined  by  'Fourier  analyzing  in  u  and  v  and  polynomial 
least-squares  fitting  in  s  the  function  values  of  Djl»Vu  and  DH*Vv 
evaluated  on  a  fixed  grid  in  the  unit  (s,u,v)  cube.  Observe  that  since 
we  already  have  a  Clebsch  representation  for  B^,  the  functions  DH«Vu  and 
D_H«Vv   can  be   evaluated  at   any    point   (s,u,v). 

The    fornula    for  DH«Vs    is   derived    from   those    for  D_H«Vu   and  DH«Vv   by 
using   the   divergence-free   condition  on  H^,    that    is, 


8  9  3 

^  DH'Vs   =    -  ;^—  DH-Vu  -  T—  DH-Vv 
8s—  3u—  3v 


The      arbitrary      function  of    u  and   v   involved    in  the    integration   of   this 
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last  equation   can   he    fixed  by  the    Fourier   analysis      of      D_H«Vs      at      some 

fixed   s. 

1      ? 
We      thus      obtain      the     Clebsch      representation     of      V   x   -yPjrB    .      A 

similar    procedure  applies    to  V   x   p^B^and    to  V   x   _B    .       If    we    dente   by     K, 

1       o 
H     and     J_  the  Clebsch   representations   of    the  vectors   V   x    P||B^   v   x   -^-pfB 

and     V      X      B^  ,       respectively,      then      the      orbit      equations      in      (s,u,v) 

coordinates   corresponding    to  equation   (4.1)    are 

4^  =  3.  Pi,  [dB   .    Vs        +        DK  •    Vsl 
d  t        m      I"-    —  —  ■' 

(4.5)  ^  =  lp,,[DB.Vu        +        DK  .    Vul 
dt        m      "  ^    —  ~ 

-^  =  i  Pii  [dB   •    Vv        +        DK  .    Vvl 
dt        m      "       —  ■" 

Similarly,    for   equation  (4.3)    they   become 

_!.  =  1  fp,.  DB   •    Vs      +      DH  •    Vs      +      ^  p?DJ   •    Vsl 
dt        m        "    —  —  Z      "    — 

(4.6)  _^  =  i  [pii  DB   •    Vu     +      DH  •    Vu     +     ^  p|?D.J   •    Vul 
dt        m   ^    "    —  —  2      "    — 

_^  =  i  [  p.,  DB   .    Vv     +      DH  •    Vv     +     ^  p^DJ   •    Vvl 
dt        m   "■     "    —  —  2      "    —  -' 

The        following        question        arises;         If         we      have      an      explicit 

representation  of      B      so      that      we      can      calculate      the      components      of 

1        0 
V    X    P||B^  ,      V      X      _     p||_B  and  V   x    B_  a  t    any   point   (s,u,v),    why  have  we    put 

these   vectors   in   Clehsch   form  before    integrating    the      orbit      equations? 

The   answer   to  this    will   exhibit   an    important    feature  of    our   method,   and 
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will  be    given  by  a   comparison  of   two    cases,    depicted    in  Figure      2a      and 
Figure  2b    of  Example   2   in  Chapter    5. 

In  Figure  2a,  we  show  the  result  of  following  a  circulating 
particle  orbit  by  integrating  equation  (4.1)  and  evaluating  the  vector 
quantity  V  x  p  B  explicitly  from  the  known  expression  for  _B  at  each 
step  of  the  integration.  Figure  2b  shows  the  same  case  when  we 
integrated  the  system  (4.5).  The  systematic  fashion  in  which  the  orbit 
expands  outward  in  Figure  2a  suggests  that  the  guiding  center 
trajectory  does  not  define  an  area-preserving  transformation  of  the 
plane  of  this  figure  into  itself,  a  fact  which  is  at  variance  with  the 
solenoidal  character  of  B_  +  V  x  p  b^  .  We,  therefore,  suspect  a 
systematic  accumulation  of  errors  during  the  integration,  and  since  a 
well-formed  drift  surface  is  clearly  indicated  in  Figure  2b,  the 
discrepancy  must  be  due  to  the  different  ways  by  which  the  right-hand 
side  of  the  differential  equations  were  computed.  The  loss  of  accuracy 
by  truncation  errors  incurred  in  computing  the  components  of  the  vector 
V  X  p  B  directly  has  destroyed  its  solenoidal  property.  These 
components  involve  the  first  and  second  derivatives  of  the 
transformation  between  physical  (x,,x9,X2)  space  and  parameter  (s,u,v) 
space,  and  so  a  major  source  of  the  errors  must  lie  in  the 
approximations  to  these  derivatives.  Moreover,  such  loss  of  accuracy 
is  inevitably  compounded  at  each  step  of  the  integration.  The  method 
we  have  adopted  to  represent  solenoidal  vectors  allows  us  to  filter  out 
the  inaccuracies  in  computing  V  x  p  B  prior  to  the  integration,  and 
permits  us  to  follow  the  orbit  effectively.  Hopefully,  this  procedure 
still   captures    the    most   essential    features    of    the    physics. 
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(b)    Scaling    che    Orbit  Equations 

The  orbit  eqijatlon  (4.1)  (or  (4.3))  appears  in  the  rationalized 
MKSQ  system  of  units.  However,  It  Is  more  convenient  to  write  this 
equation  in  terras  of  nondlmenslonal  qiiantltles  and  to  Integrate  the 
dlmensionless  form  of  the  differential  equations.  Let  us  agree  to 
express  all  dlmensionless  variables  with  primes  (')  unless  otherwise 
stated.  For  example,  B'  Is  the  dlmensionless  magnetic  field  strength. 
^1  physical  quantities  may  be  scaled  from  their  dlmensionless  values 
according  to   the    following    formulae: 

(4.7)  1  =  B^   B' 

(4.8)  ii  =  L  2i' 

(4.9)  :!L  =   Vq   v' 

where  Bq  Is  the  main  toroidal  field  strength  measured  In  teslas,  L  Is 
the  average  minor  plasma  radius  measured  In  meters,  and  Vq  Is  the 
Initial  speed  me;isured  In  meters  per  second.  With  these  scallngs,  the 
dlmensionless  time  variable  t'  and  the  dlmensionless  variable  pu  of 
Boozer    are   defined   by 


(4.10)  t      =  _L  t' 

^0 


and 

"  qBQ    B'  qBg   ^ H 

We  substitute  equations  (4.7)  through  (4.11)  Into  (4.1),  and  drop 
the  prime  notation  to  obtain  the  dimeas lonles s  orbit  equation  for 
circulating    particles   as 
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(4.12)  ^  =  .iL[B    +(Jl)V    X    p„B] 


where    the   constant 

r^  mvg 


L  qBgL 

is      the      ratio      of      the    particle    gyroradlus    to    the   minor    plasma    radius. 

For    the   guiding   center   approximation    to  be    valid,    this    ratio   must   be      a 

small   number.      Equation   (4.12)   makes   explicit   the    fact    that    the    guiding 

center   orbit    equation      is      a      small     perturbation      of      the      field      line 

equat  ion. 

\     trapped      particle      may      undergo   reflection    in    the    local   mirrors 

arising     from      tlie      helical      ripples      of      the      magnetic      field.  Such 

reflection  may    take   place    in    time    intervals    comparable    with  wq      ,    where 

Wrj  is    the  gyrof requency.      We   are,    therefore,    led    to   scale    quantities   so 

that      this      trapping      behavior  may    be  captured.      Accordingly,  we  define 

the   dimensionless   variables  _3 ' ,    x'    as  before,   and    t'   by 

(4.13)  t      =   t'/wQ 

where      Wq      =     qBQ/m     is      the      gyrof  requency      of      the        particle.  The 

nondi:Tiensional      velocity     v_' ,    energy  W,   magnetic   moment   y  '  ,    and  Boozer 
variable  p.,'   are    then  given  by 

(4.14)  V     =     L  oiQ  v' 

(4.15)  '  '  ■       W      =     m  L-^  oj^  W' 


(4.16)  p      =     mL^  0)2  y  '    /  B^ 
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(A.17)  P|l     =     L  p|| 

and    the   ratio  of    the   gyroradius    to   the  minor  plasma   radius   is 


rj  ravQ  ,  ,„ 

-r    =    ^    =    (2w')^/2 

L  qBpL 


Substituting  equations  (4.7),  (4.8),  (4.13),  (4.14)  and  (4.17) 
into  equation  (4.3),  and  then  dropping  the  prime  notation,  we  obtain 
the  dimensionless  orbit  equation  for  trapped  particles  in  the  final 
form 


(4.18)  _  =  P||B      +    V   X    l-p^B      +     y  P?  V   X    B 


The   dimensionless   form  of   the   differential  equation  (4.4)  for  the 
Boozer  variable  p.,  may  be  similarly  computed  to  be 


(4.19)     _1  =  B  .  V  i-p2   +  P|,  V  X  B  .  V  l-p2 


1      ? 
with   the   nondimensional  vector  V   —  p||    given  by 


'^,f    -(^i^jM^lVB 


Our    computer   programs    were   coded    to    integrate    the      orbit      equations      in 
their   dimensionless    forms   as    just   described. 


•40- 


5.    A  DISCUSSION   OF    SOME   EXAMPLES   AND   RESULTS 
(a)    The  Equilibria  Studied. 

The  computer  code  of  Betancourt  and  Garabedian  has  been  used  to 
produce  equilibrium  magnetic  fields  in  stellarator  geomtries.  We  have 
considered  one  case  that  was  run  with  3  =  0  and  another  with  6=4%, 
where  3  is  the  plasma  beta  defined  to  be  the  percentage  ratio  of  the 
maximum  plasma  pressure  p  to  the  maximum  magnetic  pressure  "^  12,  In 
fusion  research,  considerable  importance  is  attached  to  configurations 
that    can  support   as   high  a    value  of   3    as    possible. 

The  equilibria  that  we  considered  assumed  there  was  no  vacuum 
region  present  and,  consequently,  that  the  plasma  extended  to  the 
perfectly  conducting  outer  wall,  the  equations  of  which  are  given  by 
(cf.      equation  (2.10),    (2.11)) 


where 


ri    =  p    cos  2iTu  -  At   cos  2it(u-v) 
z.    =p    sln2iTu  +  A2sin2iT(u-v) 

p      =    [1    +  ^^^3  cos  2tt(3u-v)]"^/^ 


and 


^2    "   "'^3    "   ^''^^ 


This  corresponds  to  a  model  of  a  conventional  stellarator  that  has  a 
combination  of  an  £  =  2  and  I  =  1>  winding.  The  Fourier  coefficients  A2 
and  A3  are  a  measure  of  the  currents  in  these  windings.  There  were  12 
field  periods  to  tVie  torus,  whose  inverse  aspect  ratio  (average  minor 
plasma      radius      to      major     plasma      radius      of      the    torus)   was  0.1.      The 
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initial    pressure   profile   was  assumed    to  be    given  by  p   =   Pq(1-s    )".      The 

rotational      transform     was     determined      by  the    requirement    that    the  net 

current  I(s)    vanishes    on  each  flux   surface    [4],      For  each   case    (3=0 

and     3=4   %),    a    set   of    these   equilibrium    runs   were  made   using   grids   of 

7   X    12  X     12  points,   9  x     16  x     16  points   and   13  x    24  x    24  points      in      the 

unit      (s,u,v)      cube.         These      particular      grids      were   chosen   since    they 

permit   mesh   extrapolation:    the   corresponding  mesh   sizes   are  _     x     x 

^  6  12 

111  1^1  1  1 

,   —  X   X   and  X   X   , 

12   '    8        16        16  12       24       24 

The  equilibrium  and  stability  properties  of  this  model  £  =  2, 3 
stellarator  have  been  extensively  studied  by  Betancourt  and  Oarabedian 
[4,5,6]  and  exhibit  many  desirable  characteristics.  In  particular,  it 
appears  that  this  stellarator  may  be  stable  to  all  physically  realistic 
raagnetohyd  rodynamic  modes  and  that  equilibrium  3  limits  of  5  %  are 
achievable. 
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(b)    EXAMPLES 
EXAMPLE    1.    Magnetic    surfaces. 

In  Figure  la,  we  show  a  sequence  of  plots  of  the  flux  surfaces  at 
four  cross  sections  of  the  torus  as  they  were  produced  by  the 
equilibrium  code  for  the  case  when  0  =  0.  The  major  axis  of  the  torus 
is  to  the  left  of  each  section  and  is  not  shown.  The  flux  surfaces 
nest  around  a  single  magnetic  axis,  marked  0,  which  is  displaced  to  the 
left  from  the  minor  axis  of  the  torus.  The  outermost  flux  surface  is 
at   s  =    1.  .    :  '  •  .-.-  ,  ,. 

Figure  lb  shows  the  result  of  following  field  lines  using  our 
method  directly  on  the  magnetic  field  provided  by  the  equilibrium  code 
for  the  case  corresponding  to  Figure  la.  The  magnetic  axis  is  marked 
X.  Five  field  lines  were  followed  approximately  5  times  around  the 
torus  to  produce  these  five  surfaces.  Since  the  torus  is  made  up  of  12 
periodic  sections,  an  intersection  point  of  the  line  and  a  cross 
section  can  be  made  in  each  field  period  to  produce  a  total  of  60 
points.  The  CDC  6600  computer  at  New  York  University  took  just  over  4 
minutes    to   produce  these    five   surfaces. 

In  our  approximation  of  _B,  we  chose  either  quadratic  or  cubic 
polynomials  in  s,  and  we  found  that  employing  no  more  than  5  harmonics 
in  u  and  3  harmonics  in  v  sufficed  to  reproduce  the  flux  surfaces  of 
the   equilibrium  code. 

Figure  Ic  shows  flux  surfaces  from  the  equilibrium  code  for  the 
case  when  3=4%,  and  Figure  Id  shows  the  corresponding  flux  surfaces 
from   the   line    tracing   code. 

In  Figure  le,  we  exhibit  the  flux  surfaces  obtained  from  field 
line    tracing    using    the  extrapolated  magnetic    field    for      the      case      when 
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3=0  (cf.  Chapter  3).  The  cross  section  is  at  v  =  0,  corresponding 
to  the  first  cross  sections  of  each  of  the  previous  figures.  The 
numbers  on  each  of  the  surfaces  refer  to  the  order  in  which  the  five 
field  lines  were  integrated  by  the  code.  Observe  that  there  are  still 
closed  surfaces  enclosing  the  chain  of  islands,  but  that  their  break-up 
appears  imminent.  The  outermost  closed  surface  was  produced  by  a  line 
with  initial  position  at  s  =  0.9,  u  =  0.0,  v  =  0.0.  Field  lines 
started  at  larger  values  of  s  were  lost.  For  equilibrium  cases  such  as 
these  that  we  are  considering  where  there  is  no  vacuum  region,  we 
define  a  line  to  be  lost  if,  in  the  course  of  the  integration,  the 
value  of  s  remains  persistently  greater  tlian  some  preasslgned  value 
Sq  >  1  without  forming  a  closed  surface  --  we  chose  Sg  =  1.2  . 
Presumably,  the  last  flux  surface  assumed  by  the  equilibrium  code  at  s 
=  1  Is  near  the  separatrlx,  outside  of  which  there  are  no  more  closed 
flux  surfaces.  Field  line  tracing  using  our  method  would  provide  a 
means  of  assessing  the  validity  of  this  assumption.  However,  care  must 
be  exercised,  since  for  values  of  s  very  much  larger  than  1,  our 
approximation  to  _B  would  be  Invalid.  One  avenue  of  future  research 
would  be  to  use  the  equilibrium  code  to  compute  a  solution  for  a  case 
that  Included  the  vacuum  region.  The  free  boundary  Is  found  as  part  of 
this  solution.  We  could  apply  our  method  to  approximate  R^  In  the 
vacu^jm  region  as  well  as  In  the  plasma  region,  and  then  trace  field 
lines    In  both    regions    to   verify    the    location   of   the    free    boundary. 

In  Figure  If,  we  show  flux  surfaces  using  the  extrapolated 
magnetic  field  for  the  case  when  g  =  4  %.  IJnllke  the  case  of  3  =  0 
(Figure     le),      no    Island   structures   were   detected.      This  suggests    that. 
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for  this  conf Iguratioa,  increasing  8  may   suppress   the   occurrence   of 
large  island  regions. 


-45- 


W). 


CROSS  SECTIONS  AT  V=:   00  ,  .  25  ,  .  50  ,  .  75  ,  1/(EP«0LZj=  0.83 
MAJOR  RADIUS^  10.00  MINOR  RADIUS^  1.00 


FIGURE  la.   Flux  surfaces  from  equilibrium  code.    B  =  0 
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FIGURE  lb.   Flux  surfaces  from  field  line  tracing.    B  =  0 
B  is  unextrapolated. 
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CROSS  SEC  n  ON:  AT  V=   0  C,  .  2  ^..  ,  ,  50  ,  .  75  ,  1/(lP«0LZ)=  0.S3 
MAJOR  RADIUS=  10.00  MINOR  RADIUS^  1.00 


FIGURE  Ic .   Flux  surfaces  from  equilibrium  code.    6  =  4  %, 
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FIGURE  Id.   Flux  surfaces  from  field  line  tracing.    6  =  -i  '= 
B^  is  unextrapolat  ed  . 
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FIGURE  le.   Flux  surfaces  from  field  line  tracing.    3=0. 
B^  is  extrapolated. 
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FIGURE  If.   Flux  surfaces  from  field  line  tracing.    g  =  A  % 
B^  is  extrapolated. 
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EXAMPLE  2.  The  effect  of  not  preserving  the  divergence-free  condition 
on    solenoidal   vectors. 

In  this  example  we  present  evidence  to  support  our  belief  that 
preserving  the  divergence-free  condition  of  solenoidal  vectors  is 
necessary  for  obtaining  flux  surfaces,  if  these  exist,  by  numerical 
integration    techniques    (cf.      Chapter    4). 

We  first  followed  the  guiding  center  orbit  of  a  circulating 
particle  by  integrating  equation  (4.12)  and  evaluating  the  right  member 
of  that  equation  explicitly  from  the  known  representation  for  B,  The 
result  is  shown  in  Figure  2a.  We  then  followed  the  same  orbit  by  first 
expressing  the  term  V  x  p^-B  of  that  equation  by  our  solenoidal  Clebsch 
representation.  The  result  is  shown  in  Figure  2b.  As  discussed  in 
Chapter  4,  the  systematic  expansion  of  the  orbit  of  Figure  2a 
manifestly  demonstrates  that  the  transformation  of  the  plane  of  this 
figure  into  itself  defined  by  the  orbit  is  no  longer  area-preserving. 
In  each  case  the  orbit  was  followed  approximately  six  turns  around  the 
torus.  '«fe  see  the  effect  of  truncation  errors  in  the  integration 
scheme  after  a  much  longer  time:  Figure  2c  shows  this  effect  on  the 
orbit  depicted  in  Figure  2b  when  we  followed  it  100  times  around  the 
torus.  We  are,  tlierefore,  led  to  the  conclusion  that  the  consequence 
of  not  preserving  the  divergence-free  condition  is  far  more  serious 
than  that  of  truncation  errors  due  to  the  numerical  scheme.  The 
integration  by  the  orbit  of  Figure  2c  used  a  step-size  At  =  0.01. 
Figure  2d    shows    the   same    orbit   as    that   of    Figure  2c   when  At   =    0.005. 

The  orbit  calculations  using  our  method  takes  roughly  twice  as 
much  time  as  a  comparable  field  line  calculation.  The  orbit  depicted 
in  Figure  2b    took    2    1/2  minutes    to  execute  on    the   CDC    6600  computer. 
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FIGUKE  2a.   Effect  of  not  preserving 
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FIGURE   2b.         V    •    K   =    0 


-5-4- 


FIGURE  2c.   Integration  step  size   At  =  0.01 
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FIGURE  2d.   Integration  step-size   At  =  0.005 
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EXAMPLK    3.    Drift   surfaces   and    the    dependence    on    the    gyroradius. 

In  this  example,  we  studied  the  dependence  of  the  drift  surfaces 
on  the  gyroradius  for  the  Jl  =  2, 3  stellarator  which  was  described  in 
Section  a  and  whose  magnetic  flux  surfaces  were  investigated  in  EXAMPLE 
1.  We  followed  the  ,q;uiding  center  orbits  of  circulating  particles  using 
the   extrapolated   magnetic    fields. 

\s  will  be  recalled  for  the  case  with  3=0,  depicted  in  Figure 
le,  the  flux  surfaces  had  island  structvires.  When  the  gyroradius  was 
1.6    %  of    the   minor    plasma    radius,    that    is,  —  =  0.016  in     the      notation 

Li 

of  Chapter  4,  the  drift  surfaces  showed  no  islands  (Figure  3a). 
%wever,  the  Islands  reemerged  when  the  gyroradius  was  increased  to  4  % 
of  the  minor  plasma  radius,  as  may  be  seen  in  Figure  3b.  The  outermost 
drift  surface  is  quite  indented.  When  —  =  0.06,  the  outer  portions 
were  populated  by  islands  only,  and  there  did  not  seem  to  be  any  drift 
surface   enclosing    them. 

For  the  case  wi th  3  =  4  %,  we  found  no  breakup  of  the  drift 
surfaces  for  various  values  of  the  gyroradius.  Figures  3c  and  3d  show 
some  drift  surfaces  when  —L  =  0.016  and  0.04,  respectively.  This  leads 
us  to  believe  that  if  the  flux  surfaces  show  no  sign  of  breaking  up 
into  islands  (Figure  If),  then  the  drift  surfaces  also  tend  to  be 
nested  for  higher  values  of  the  gyroradius  than  if  the  flux  surfaces 
had  island  structures.  \r  with  the  flvix  surfaces  for  this  £  =  2,3 
stellarator,  higher  values  of  3  tends  to  suppress  the  occurrence  of 
islands    in  the    drift    surfaces. 

In  Figvire  3e,  we  show  the  deviation  of  drift  surfaces  from  a 
magnetic  surface  as  the  gyroradius  is  varied.  The  magnetic  field  used 
was   uneKt  rapolated .       In    this    plot,    the    outermost   surface    (labeled    2)    is 
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supposed  to  be  the  plasma/ vacuum  interface  and  is  included  merely  for 
reference.  Three  other  surfaces  are  shown:  one  is  a  flux  surface,  and 
the  other  two  are  drift  surfaces  corresponding  to  gyroradii  that  are 
3  %  and  6  %  of  the  minor  plasma  radius.  Each  surface  was  produced  by 
starting  the  trajectory  at  the  same  point  in  space,  where  the  poloidal 
angle  was  —  (u  =  0.125).  Proceeding  from  the  magnetic  axis  outward  to 
the  free  boundary  along  the  ray  u  =  0,  we  encounter  the  magnetic  flux 
surface  (labeled  1)  first,  then  the  drift  surface  corresponding  to 
gyroradius  3  %,  and  then  the  drift  surface  correponding  to  gyroradius 
6  %  (labeled  4).  In  relation  to  the  flux  surfaces,  the  drift  surfaces 
tend    to   be    shift    outwards    from   the   major  axis   of    the    torus. 
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FIGURE  3b.      Drift    surfaces.  6=0.      -  0.0-4 
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FIGURE  3c.  Drift   surfaces,    B  =  A  "/: .  ~  =   0,016 
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FIGURE  3d.   Drift  surfaces.    6=4%.  -^  =  O.OA 
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FIGURE  3e.   Deviation  of  drift  surface  from  magnetic  surface, 
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EXA-tPLR    4.    Trapped    particles. 

In  this  final  example,  we  demonstrate  the  ability  of  our  method  to 
follow  the  orbits  of  trapped  particles  using  the  £  =  2,3  stellarator 
configuration   described   in   Section   a    of    this    chapter. 

Figure  4a  shows  a  plot  of  the  orbit  of  a  helically  trapped 
particle  in  the  poloidal  (r,z)  plane  of  the  torus.  In  the  literature 
[16],  such  an  orbit  is  called  a  helical  banana.  The  particle  bounces 
between  the  tips  of  its  baaana  orbit,  and  the  banana  as  a  whole 
precesses  clockwise  around  the  minor  axis  of  the  torus.  The  tips  of 
the  banana  are  reflection  points  of  the  orbit  where  the  parallel 
component,  v.,  ,  of  the  particle's  velocity  vanishes  (cf.  equation 
(4.2)).  The  time  interval  betvjeen  successive  reflections  is  called  the 
bounce    time. 

Figure  4b  shows  the  intersection  of  the  particle's  trajectory  with 
a  fixed  toroidal  plane,  v  =  0.7.  Here  the  banana  shape  is  more  evident 
and  is  not  unlike  a  magnetic  or  drift  surface  island  of  some  of  the 
previous  examples.  The  time  taken  for  the  particle  to  bounce  across 
the    widest    part    of    this    banana    is    a    good    measure    of    the    bounce    time. 

We  integrated  the  dime nsionles s  guiding  center  equations  (4.8)  and 
(4.9)  to  obtain  this  orbit.  *is  we  mentioned  in  Chapter  4,  for  a  given 
magnetic  field  configuration,  there  is  a  t wo-par.^me ter  family  of  orbits 
depending  on  the  energy  W  of  the  particle  and  on  its  magnetic  moment  y . 
Since  the  differential  equations  are  dime ns ionless ,  some  guide  as  to 
what  magni  bides  to  use  for  the  dime  as  ionless  energy  and  magnetic  moment 
is  neeiied .  To  this  end,  we  applied  the  laws  defined  by  the  equations 
(4.7),  (4.8)  and  (4.13)  through  (4.17)  to  scale  the  physical  model  of 
following      the       orbit      of       a      30   keV      proton      in      a      12      period    £    =    2,  3 
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stellaraCor  with    inverse  aspect    ratio  of   O.l   and  average  minor      plasma 

radius      of      50     cm;       the     main      toroidal      magnetic    field  strength  B  was 

assumed    to  be    1.6   teslas.      The  gyrof  requency ,   Wq   =  qB/m,  is   then    1.5     x 

Q 

10  Hz,  and  tlie  dime ns ionles s  energy  (cf.  equation  (4.15))  and 
magnetic  moment  (cf.  equation  (4.16))  are  0.5  x  lO"-*  and  0.5016  x 
10  -' ,  respectively.  It  is  precisely  because  these  last  two  values  are 
so  nearly  equal  ttiat  the  trapping  depicted  in  Figure  4a  occurs  (cf. 
equation  (4.2)).  Indeed,  when  we  lowered  the  tragnetic  moment  to  0.37  x 
10  -" ,  the  particle  was  not  trapped  and  began  to  circulate  as  shown  in 
Figure  4c.  We  followed  the  orbit  of  Figure  4a  over  1000  gyroperiods 
(=  1000  Wq  ) ,  and  since  60  reflections  are  Indicated  In  that  Figure, 
tlie    bounce    time    Is    approximately    16.7  gyroperiods. 
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FIGURE  4a.   Helically  trapped  particle  Cr,z)  plane.    W  =  0.5  y   lo" 
y  =  0.5016   X  10"'^. 
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FIGURE  4b.   Helically  trapped  particle  at  v  =  0.7 
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FIGURE  Uc.      No   trapping.      W  =   0.5    x   lO"'^;      u   =    0.37    x   lo""^ 
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(c)    Discussion 

Collisions  between  charged  particles  in  a  plasma  constitute  one  of 
the  principal  causes  of  the  diffusion  of  these  particles  and  the 
subsequent  transport  of  energy  across  field  lies.  It  is,  therefore, 
insufficient  to  study  plasma  confinement  solely  on  the  basis  of  the 
geometry  of  the  flux  surfaces,  drift  surfaces  or  orbits  of  isolated 
trapped  particles  not  subject  to  collisions  [14],  Nevertheless,  such 
calculations  are  still  important  even  in  the  context  of  a  sophisticated 
theory  of  diffusion,  since  they  reveal  the  possible  location  of  a 
particle  between  collisions.  Suppose  a  particle  is  circulating  in  a 
region  where  the  drift  surfaces  are  nested,  then  a  collision  between 
this  particle  and  another  will  cause  it  to  scatter  a  radial  distance, 
Ar,  onto  another  drift  surface.  Since  the  surfaces  are  assumed  to  be 
nested,  the  effective  radial  excursion  of  the  particle  after  two 
collisions  can  only  be  as  much  as  2Ar.  However,  if  the  particle  is 
circulating  in  a  region  where  the  drift  surfaces  contain  substantial 
islands  of  width  w,  very  much  greater  tlian  Ar,  then  the  effective 
radial  excursion  of  the  particle  after  2  collisions  can  be  as  great  as 
w  +  Ar.  In  general,  one  may,  therefore,  say  that  the  presence  of 
islands  in  the  magnetic  surface  or  drift  surface  structure  will  lead  to 
enhanced  particle    diffusion  and   energy   transport. 

The  situation  with  trapped  particles  is  similar.  Consider  the 
helically  trapped  particle  orbit  depicted  in  Figures  4a  and  4b.  As  we 
have  mentioned,  the  particle  traverses  the  width  w  of  the  Island  shown 
in  Figure  4b  over  one  bounce  period.  If  the  frequency  of  collisions 
among  particles  is  so  low  tliat  the  time  interval  between  collisions  is 
comparable      to      the    bounce    time,    then    the   effective    radial   excursion   of 
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the  particle  may  be  increased  by  the  amount  w.  Such  considerations  are 
crucial  for  an  adequate  evaluation  of  diffusion  and  transport 
especially  in  stellarator  devices  operating  in  low  collisionali ty 
regimes.  Although  a  large  body  of  theoretical  knowledge  already  exists 
in  the  theory  of  neoclassical  diffusion  [16],  many  important  questions 
remain  largely  unanswered.  One  such  is  certainly  the  effect  of  3  on 
diffusion.  Tt  is  our  hope  that  the  method  described  in  this  work,  when 
used  in  conjunction  with  the  Be tancourt-Garabedian  code,  may  contribute 
to   the    resolution   of   these    issues. 
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6.    THE    COMPUTER    PROGRAMS    COEFF   AND   LINORB 

We  have  included  tlie  listings  of  two  FORTRAN  programs,  COEFF  and 
LIMORB,  in  Sectin  b  of  this  chapter.  The  two  programs  implement  our 
method  of  representing  the  nvignetic  field  and  of  following  field  lines 
or  guiding  center  orbits.  Both  codes  are  internally  documented,  and  so 
we  shall  only  describe  their  operation  as  well  as  some  pertinent 
features   in    this    chapter.  ...,.       ,  , 

Since  our  method  is  predicated  on  being  in  possession  of  a 
numerical  equilibrium  solution  from  the  Be tancourt-Garabedian  code, 
BETA,  it  should  be  clear  tlvTt  before  either  COEFF  or  LINORB  can  be 
used,      an     adequate      equilibrium     has      been      computed      by      BETA.  The 

documentation  of    BETA  appears    in    [3]. 

(a)   Operation   of   the   Programs 

The  program  COEFF  computes  the  series  represetations  of  the 
components  of  the  vectors  _B  and  K  =  V  x  p..  B  as  described  in  Chapters  3 
and  4,  respectively.  In  addition,  COEFF  also  computes  the  series 
representations  of  the  functions  rQ(v),  Zq(v)  and  R(s,u,v)  which  appear 
in  the  mapping  between  the  unit  (s,u,v)  cube  and  physical  (r,z,e)  space 
(cf.      equations   (2.9),    (2.10),    (2.11),    (3.11),    (3.12),    (3.13)). 

Ihe  program  LINORB  follows  field  lines  and  guiding  center  orbits 
of  circulating  particles  by  integrating  the  system  of  equations  (4.12). 
To    follow    field  lies    the    gyroradius    is    merely  set    to   zero. 

Both  progrims  liave  been  designed  to  allow  the  user  to  submit  runs 
as  BATCH  jobs  from  a  remote  terminal.  Our  codes  have  been  run  on  the 
CDC  6600  computer  at  the  Courant  -la thema tics  and  Computing  Laboratory 
of  New  York  University.      They    make   use   of    two    library    routines    that   may 
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be  specific  to  this  facility  only.  Consequently,  some  alterations  may 
be  necessary  if  our  programs  are  used  on  other  computer  systems.  We 
refer  to  the  subroutine  LEQ  and  to  the  Calcomp  plotter  package,  SCOOP. 
LEQ  is  a  system  library  routine  used  to  solve  a  linear  system  of 
algebraic  equations;  it  is  called  by  the  subroutines  SIMT,  SINTRO  and 
CSI^T  of  the  program  COEFF.  Since  most  computing  facilitines  support 
software  packages  that  have  similar  subroutines  to  LEQ,  the  replacement 
of  the  CALL  statements  invoking  LEQ  should  present  no  difficulties. 
Indeed,  in  the  past,  we  have  successfully  used  the  subroutines  LEQTIF 
in   the  IMSL  package,    which   is    widely   supported. 

The  library  package  SCOOP  is  invoked  by  calls  to  the  subroutines 
PLOTSBL,  PLOT,  SYMBOL,  NUMBER  and  FRAME  which  appear  in  the  subroutine 
GRAFIX  and  0RBPLT3  of  our  code  LIMORB.  The  subroutine  PLOTSBL 
initializes  the  SCOOP  plotting  package.  Calls  to  subroutine  FRAME  set 
the  pen  on  a  new  page,  as  required.  Actual  pen  motions  within  each 
page  are  controlled  by  calls  to  the  subroutines  PLOT  and  SYMBOL. 
Subroutine  NUMBER  is  invoked  merely  for  labeling  purposes  and  is  not 
essential.      The   last    CAIJ.   to  subroutine    PLOT  appearing    in  GRAFIX, 

CALL   PLOT(0. 0,0. 0,999) 
is    required    for   normal    termination    of    the    SCOOP    package. 

At  the  concl\ision  of  one  successful  run  of  the  equilibrium  code, 
BETA,  a  data  file,  TAPE  7,  is  created  and  must  be  saved  since  it  is 
used  as  input  to  the  program  COEFF.  TAPE  7  contains  the  computed 
values  of  the  functions  Tq^zq,  R  and  z,  described  in  Chapter  2.  These 
data  are  written  onto  TAPE  7  just  before  the  termination  of  one  run  of 
BETA.         The   format    of    the    data    on   T.APE    7   may    be   deduced    by    locating    the 
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READ(7,«)   statements    in   the   following  subroutines    of    the    program  COEFF: 
READl,     SUBAX,    SUBRO,    SUBPSI. 

In   addition    to  TAPE    7,    the    program  COEFF  requires    that   a   number   of 
parameters     be      initialized        properly        before        compilation.  These 

parameters  are  defined  in  FORTRAN  DATA  statements  in  the  subroutine 
UNIV.  The  definitions  of  these  parameters  are  independent  of  the  data 
on  TAPE  7.  This  permits  the  user  to  vary  their  values  using  the  same 
TAPE  7  from  the  equilibrium  code.  The  meanings  of  these  parameters  are 
described  in  subroutine  UNIV.  For  the  most  part,  they  refer  to  the 
degree  of  the  approximating  polynomials  in  s  (IMAXRO,  IMAXPSI,  IDEG), 
the  number  of  harmonics  in  u  (JMXRO,  JMXPSI,  .IMAX)  and  the  number  of 
harmonics  in  v  (KMXRO,  KMXPSI,  KMAX,  KAXIS)  used  in  the  representations 
of  B,  K  and  the  transformation  functions.  After  the  initialization  of 
these  parameters,  COEFF  may  be  compiled  and  executed.  Fourier  analysis 
of  the  data  from  the  equilibrium  code  is  carried  out  by  a  Fast  Fourier 
Transform  subroutine,  FFORM.  Approximations  in  the  variable  s  employ 
least    squares    fitting  of    data    to  a    polynomial 
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This  takes  place  in  subrouties  SINTRO,  SINT  and  CSIOT.  A  data  file, 
TAPE  83,  is  created  at  the  end  of  a  run  of  COEFF  and  is  used  as  input 
to  the  code  LINORB.  TAPE  83  contains  the  coefficients  of  the 
aforementioned  series  of  the  vector  components  and  transformation 
functions.  The  format  of  the  data  on  TAPE  83  may  be  found  by  a  perusal 
of    the    subroutine   READl    of    the    code    LINORB. 

As     with     COEFF,      the      program        LINORB        requires,        for        proper 


-73- 
i nitiallza tion ,  Che  presence  of  a  data  file  (TAPE  83),  and  the 
definitions  of  various  parameters,  which  are  independent  of  the  data  on 
T.\PE  83.  These  definitions  are  made  in  DATA  statements  in  the  main 
routine  of  LII^RB.  Their  meanings  are  described  there.  The  parameters 
Include  the  number  of  lines/orbits  to  be  followed  (NORB),  the  number  of 
times  around  the  torus  each  line  or  orbit  is  to  be  followed  (NVMAX), 
their  starting  positions  ( STV) ,  and  the  integration  step-sizes  for  each 
trajectory  (DTO).  \s  listed  below,  LINORB  will  integrate  up  to  5 
trajectories.  The  solution  of  each  integration  Is  stored  on  a  separate 
tape  which  may  be  saved  for  future  use:  trajectory  £  Is  stored  on 
TAPE  £.  The  integration  procedure  uses  a  fourth-order  Adams-Moulton 
nultlstep  scheme  with  a  fourth  order  Runge-Kutta  singlestep  starter. 
The  program  also  allows  for  halving  or  doubling  the  integration  step 
size  depending  on  whether  the  local  truncation  error  Is  greater  or  less 
than  preasslgned  values,  respectively.  The  Integration  of  a  trajectory 
will  stop  If  (1)  the  trajectory  has  traversed  the  torus  the  designated 
number  of  times,  (2)  the  trajectory  goes  beyond  s  =  1.2,  or  (3)  the 
mesh-halving  subroutine,  HALVE,  has  been  invoked  ten  consecutive  times. 
Diagnostic  messages  are  printed  at  the  conclusion  of  each  run  to 
Indicate  how   each    Integration   was   terminated. 

Finally,  ttie  plotting  of  the  solutions  takes  place  after  all  the 
trajectories  have  been  Integrated.  Planar  cross  section  puncture 
plots,  such  as  Figiire  le  of  Chapter  5,  are  produced  by  the  subroutine 
OR3PLT3.  The  user  chooses  at  which  plane,  v  =  constant,  within  a  field 
period  the  plots  are  desired  by  appropriately  Initializing  the  variable 
VPL  before  running  tlie  code.  In  general,  since  we  shall  not  have  the 
value    of    the    solution   on    precisely   the    plane      at      which      we      desire      to 
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plot,      tVie    value  of   the    solution  at    the   cross   sectional    plane   has    to  be 
approximated.      We  do    this    by  a    linear    interpolation   between   two      values 
of    the   computed    solution    that    are   nearest   to   that    plane. 

The  operation  of  our  codes  COEFF  and  LINORB  is  quite  simple.  The 
number  of  parameters  at  our  disposal  to  vary  is  small.  Both  codes  have 
been  written  largely  in  block  form  to  facilitate  alterations  where 
necessary.  It  is  hoped  that  these  programs  operated  in  conjunction 
with    the  parent    equilibrium  and    stability    code    will   prove   useful. 
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N  THE  S 
INE  HAL 
bOJNJ 
ERROR 
BY  A  C 
1^2^3,^ 
RUNC AT  I 
GTH  SCA 
/IN  SU3 
THEN  XL 
IVE  NUM 
RETURN 


RBITS 
f  NOR 
INE  L 
US.  N 
ATER 
EP  SI 
BIT. 
TER  T 
U  I  fM  T 
LABE 
EATER 
)  THE 
T  WHI 
OLERA 
OR.  I 
TEP  S 
VE. 
ON  TH 
IS  LE 
ALL  T 
)  WE 
ON  ER 
LE  US 
ROUTI 
ENGTH 
BEP;  I 
RLARM 


TO  B 
B  CAN 

IS  T 
VMAX 
THAN 
ZE  TO 
DTO  I 
HAN  0 
ES  OF 
LLED 

THAN 

V  CO 
CH  PO 
TED  V 
F  THI 
IZE  I 


E  fOLLO 

TAKE  G 
0  BE  FO 
IS  A  VE 
OR  EJUA 

BE  USE 
S  A  VEC 
R  EQUAL 

THE  IN 
L  .THE 

OR  ECU 
QRCINAT 
INTS  WI 
ALLE  OF 
S  ERROR 
S  HALVE 


E  LOCAL  TRUN 
SS  THAN  SEP, 
0  THE  SUBRUU 
IGHTS  USED  I 
ROP. 

ED  TG  CCMPUT 
NE  COLARM.IF 
(L)  MUST  BE 
N  THIS  CASE 
=0.  OTHER WIS 


WED. 

N  VALUES  FROM  1  TO  5. 

LLOWED 

CTDR  ARRAY  WITH 

L  TO  NORB. 

D  IN  INTEGRATING 

TOR  ARRAY  WITH 

TO  NORB. 
ITIAL  POSITION  OF 
COLUMN  DIMENSION  OF 
AL  TO  NORB. 
ES  WITHIN  A 
LL  BE  PLOTTED. 

THE  LOCAL 

IS  GREATER  THAN 
0  BY  A  CALL  TO 

CATION  ERROR. 
THE  STEP  SIZE  IS 
TINE  DOUBLE. 
N  COMPUTING  THE 

E  THE  GYRORADIUS* 

L  IS  A  FIELD 
SET  TO  ANY 
SUBROUTINE  COLARM 

E  »  R  L  A  R  M  IS 
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C 
C 

c 
c 
c 
c 
c 
c 
c 


(IDBi JDB>KD3) 


(  IDK, JDK^KDK 
(  IDk,  JDR>KDf^. 


CDMfl 
COMM 

+ 
+ 

+  ss 
+  cc 
+     cc 

C  G  ;<  •- 
CDMM 
COMM 

+ 
COM'^ 
COM,^ 
COMM 
COKM 
COMM 
CO;iM 
COMM 
DIME 

+    S 
DATA 

+ 
DATA 

+ 

+ 

+ 

+ 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
PI2- 
XPE»^ 
DONE 
REA 
CALL 
PRI 
WRIT 
WRIT 


ON/VAPl/P 
QN/VARIO/ 


DN/VARll/ 
1K(6^10*1 
2K(6,10,1 
3K(6,10>1 

■"  \  /  V  A  f  i  5  / 
DN/VAR131 
OM/VARl^/ 

0N/\/Akl5/ 
QN/VAR17/ 
uN/yAP171 
CN/VAkl8/ 
0N/VAR19/ 
QN/VAR20/ 
0N'/VAR21/ 
MS  I  ON  X(^ 
TV(3f 5) >H 
NQP6,NVM 


STV 


VPL 

GEP 

WGT 

XLE 

IDB 

IDK 

IDR 

2.0* 

=  1.0 

'-99 

D  IN 

REA 

NT  I 

E(6, 

E<6i 


/O.bf 
0.7» 

o,e» 

0.5/ 
0.33 
/0.02 
/SEP/ 
/l.O, 
NGTH/ 
,  JOB* 
»  JDK/ 
/  JDR/ 
ACQS( 

9.0 

DATA 
01 

NPUT 
597) 
600)Q 


INVER 
IS    A 
OR     EQ 

TH 

TH 
)     '     TH 

TH 
)     =     TH 

TH 
I2/EP/ 
CC1(6/ 
CC2(6/ 
CC3(6/ 
CC1K(6 
0), 

0)/CS2 
0),CS3 

r  ■■,  J" 

/IMK/J 
RCC  (6/ 
IRO/JR 
CSU2/S 
RLARM 
/XLENG 
P( 15)/ 
R  A  C  {  2  6 
DELO/D 
V  P  L  (  ^  ) 
)/Y('^) 
(  '.  )  /  H  X 
AX/DTO 
.01/  . 
0.0/ 0, 
0.0/C. 
0.0/ 0. 
0.  0/ 0. 
33/0.0 
/0.27/ 
.lE-2/ 
1.0/1. 
500.0/ 
KC6/6/ 
KDti/b, 
KDR/6/ 
-1.0) 


SELY  P 
VECTOR 
UAL  TO 
E  DIME 
E  CliEF 
E  DIME 
E  CJEF 
E  DIME 
E  CUEF 
QLZ/ AL 
10/10) 
10/10) 
10/10) 
/lO/lO 


RQPORTI 
ARRAY 
NORB. 
NSICNS 
FICIENT 
NSIONS 
FICIEfiT 
NSIONS 
FICIENT 
F 

/CS1(6/ 
/  C  S  2  (  6  / 
/CS3(6/ 
)/CSlK( 


ONAL     TO    XLENGTH(L).     XLENGTH 
WITH    DIMENSION    GREATER    THAN 

OF     THE     ARRAYS    THAT     CONTAIN 
S    OF     DEGRADS/DBGRADU/DBGRADV. 
OF    THE    ARRAYS    THAT    C0r4TAIN 
S    OF     DKGRADS/OKGRADU/DKGRAOV. 
QF    THE     ARRAYS    THAT    CONTAIN 
S    OF    RG    IN    THE    MAPPIImG. 

10/10)/SCl(6/10/10)/SSl(o/10/10)/ 
10/10)/SC2(6/10/10),SS2(b/10/iO)/ 
10/10)/SC3{6/ 10/10) /SS3(b/lO/ 10) 
6/10/lC)/SClK(6/10/10)/ 


K{6/1C/10)/SC2K(6/1C/10)/5S2K(6/10/10)/ 
K(6/10.10)/SC3K(6,1C/1C)/SS3K((;/10/10) 

/  K  [;  ^ 

K/KDK 
10/10)/RSC(6/10/10)/RSS(6/10/10)/ 

/KDR 


J  t  \r.b/  I  b'-i/  Jlif 

MK/KMK/  IDK/ JD» 
1C/10)/KCS (6/ 

0/KkD/ IDR/ jdr, 
INU/CSV2/SINV 


TH(5 

Q(  15 

)/RA 

ELI/ 

/  NOR 

/SU'^ 

(^/o 

/5/3 

01/  . 

0/ 

0/ 

0/ 

0/ 

/CO 

0.52 

.lE- 

0/0. 

5C0. 

10/  1 

10/1 

10/1 


) 

)/PP(15>/ 
S(26)/ZAC 
DtL2/DEL3 
B/XPER 
(^)/PSUM( 
)/ XX(^/6) 
0/30/30/3 
01/  .Ol/.O 


Q  Q  {  1  5  )  /  R  R  (  6  ) 

(26)/ZAS{26)/KRAX/KZAX 

/DEL10/DEL20/DEL30/DEL22/DtL33 

A)/CSUK(^)/DEL( A )/WGT(4)/ 

/YP(^)/YC('i)/DT0(5)/NVNAX(5) 

0/30/ 

1/ 


,0,11/ 
5/ 

0/ 

0/500. 0/^00. 0/^CO.O/ 

0/ 

0/ 

0/ 


AND     INITIAL    PARAMETERS 
LZ/tP/ALF/DEL0/DELl/DEL2/CEL3/DEL10/DEL20/DEL30/ 
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602 


WRITE 

WRITE 

597  FORMA 

600  FCR-IA 

+ 

+ 
■f 
+ 

601  FORMA 
+ 

+ 

+ 

603  FORiA 

DO  90 

R  E  W  I  N 

C   COMPU 

CALL 
C   I N  I  T  I 
WRITE 
FORMA 

+  *RLAf^ 

+10X, * 
VMAX  = 
OTMAX 
D  T  M  I  N 
NHALV 
NOOUB 
NSTEP 
DO  1 
XX(  II 
XX(  II 
HX(  II 
.  HX{  II 
DT   «= 
X  (1) 
X(2) 
X(3) 
X  (  ^) 

GeTAi 

DC  10 

DO  11 

s  u  r  {  I 

CALL 
DO  12 
HX(I  I 
D£L(  I 
SUM(  I 
CALL 


(6, 

(6, 

T(l 

T(  / 

5 

5 

5 

5 

5 

T(  / 

5 

5 

T(  / 

L  = 
D  1 
TE 
COrt 
ALI 
(6, 
T{  / 
M«* 
S  =  * 
QLZ 
=  2. 
'0. 
=  0 
»0 
»1 
11  = 
^6) 
>5) 
f  6) 
>5) 
DTJ 


11 


12 


\    3 

1  = 

I  I 

I)  =• 

EVA 

II 

»5- 

I)- 

I  )« 

EVA 


DO  13  .11  = 


D 

601)  I 
6C3)G 
H1/2X 
/SX>* 
X,*  I'i 
X,*AL 
X>*DE 
X,*DE 
X,*DE 
/5X,* 

x,*n 

X,  *i<R 
X>*IM 

/bx,* 

l^NOR 
0 

RLARH 
LARM( 
ZE  EA 

602)  L 
/5X,* 
>E20. 
»F10. 
*NVr,A 
0+DTC 
01«DT 


1*^ 

•DONE 
"DOME 
=  C  G  iN  E 
=  DOHE 
(L  )/4 
=  STV{ 
=  STV( 
=  STV( 
»0.0 

STAk. 
1/  3 

0  .0 
L{X(1 
=  1,^ 

1  )«H( 
DT*H{ 
SUM  (I 
L  (X(l 

X(^) 
I*'. 


EL22*DEL33 

MBfJMB^KMtiWROiJPOfKRO^KRAXfKZAX^IMK/JMKiKMK 

EP/SEP 

»  ♦  f-  E  S  L'  L  T  5  *  ) 

CLZ=*>Flp.5/ 

VERSE  .•'AJJR  RADIUS  OF  TORU  S  »«  /  F  10.  5  / 

f=*, FIO. 3/ 

LO  «*»F10.  5/2X>*0tLl  =  *  *  Fl  0  .  f)  >  2  X^  ♦  DE  L2  =**F10.5/ 

L3  =*,F10.5^2X,*DfL10  =  *>F10.3/2X,*DEL2  0  =  *,F10.5>/ 

L3  0=*>F10.ti>2X»*DEL22  =  *>F10.5»2X/*DEL33  =  *^F10.5) 

1MB  =**I5>2X»*JMB  =**I5^2X**KMB  =*»I5/ 

G  ='^>I5i2X^«JPa  ==f,I^^2X>«KR0  »->,I5/ 

AX=*,  I5^2X,*KZAX  =  *,  lb/ 

K  «*»I5/2X/«JMK  «*>  I5*2X**KMK  »*^I5) 

&EP   =*>E20,10/5X»*SEP   =*^E20.10) 

8 


RLARMi  L  ) 

CH  Of-BIT.  .  . 

,  XLENGTh(L)»RLAF<H>STV(l>L)>STV(2>L),STV(3^L)^DT0(L) 

OP  BIT  *jI3/7X/*XLENGTH  =  *>F12.5*2Xf 

1G/7X>*I.NITIAL  POSITION:*/ 

5>2X,*U=*iF10.P^2X,*V=*>F10.5>2X,*0T=*^F10.5) 

X(  L) 

(L  ) 

0(L  ) 


.0 

1;L  )<'*ALF 

2>L  ) 

TING  VALUES  BY  Tht  RUNGE-KUTTA  METHOD. 


)*X(2)*X(3)^X(^)/H) 

II) 

II) 

I  )  +  DEL (  II  ) 

)+0.5*DEL(l)>X(2)  +  0.5*DEL(2)>X(3)+0.t)*DEL(3), 

+  0.t.*DEL(^)*H) 
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13 


l*. 


15 


10 


16 


162 

C 

C 

C 

C 
161 

17 

18 

19 

20 


21 


22 


DEL(  I 
SUM(  I 
CALL 

DO  1^ 
DEL(  I 
SUM(  I 
CALL 

DO  15 
DtL(  1 
SUM(  I 
Y(  II  ) 

iF{r( 

DO  10 
XX(  J> 
X(  J)» 

CONTl 

CALL 

DO  16 

HX(  II 

XXdl 

DT  =  DT 

NSTEP 

DO  It 

WRITE 

CDNTI 

kUNGE 

THE  A 

THE  P 

CALCU 

DO  17 

PSUM( 

CSUM( 

DO  lb 

PSUMC 

DO  19 

PSUM( 

CSUM  { 

DO  20 

PSUM( 

CSUM( 

CO  21 

PSUM( 

CSUrtC 

YP(  II 

CALL 

DO  22 

CSU.1( 

YCdl 

IF(YC 

CALL 


I)='DT*H( 

I  )-5UM(l 

EVAL{X(1 

X(^) 

II'l*^ 

I )-DT»H( 

I  )  "SUM (I 

EVAL(X(1 

X(^) 

II«1*A 
I  )»DT*H( 

1  )  =  S  U  M  (  I 
«X(  II  )+S 
1)  .LT.O. 

J»l,4 
5-1  )  =X( J 
Y(  J  ) 
NUE 
EVAL(  Yd 

II=1>4 
>1) «H(  II 
d)«Y(  II 
0(  L) 
=  5 

2  I'lf'f 
(10)  XX( 
NUE 

-KUTTA  I 
DA'^S-MOU 
-  AND  C- 
LATING  T 

II  )'=0.0 
II  )  =C.O 

I  I  =  l>^ 
II  )=PSUM 

II  =  l>'t 
ID^PSU-i 
II  )«CSUM 

11=1/^ 
II) «PSU1 
II  )  «=CSU^l 

II-l^'. 
II  )  =  P  S  LI .", 
II  )-CSUM 
)»XX( 11^ 
EVAL( YP( 

II«1,<. 
II  )  «=CSUM 
)«XX(  II, 
d)  .LT.O 
&VAL(YC( 


II) 

I  )  +  2.0*DEL(II ) 

)+O.5*DEL(l)»X(2)  +  0.5<'0EL(2),X(3)  +  0.5*DEL(3), 
+0.5*DEL{^)*H) 

II  ) 

I  )+2.0*DEL{ II ) 

)+DEL(lJ^X(2)+DEL(2),X{3)+DEL(3)> 
+DEL 14)*H) 

II  ) 

I)+UEL(  II) 
UM(  II  )/fc.O 

O.OR.Yd)  .GT.  1.2)  GO  TO  89 


)>Y(2)/ Y(3), Y(4),H) 


l»5-I)^XX(2i5-I)>XX(3>l?-I),XX(^,5-I) 

S  NOW  DONE.  THE  NEXT  STEP  IS  OBTAINED  BY  USING 
LTCN  PREDICTOR-CORRECTOR  MULTISTfcP  METHOD. 

PREFIXES  REFER  TO  QUANTldES  INVOLVED  IN 
HE  PREDICTED  AND  CORRECTED  VALUES^RESPECdVELY. 


(  II  )-9.C*HX( II>^ ) 

(  II  )+37.0*HX(  11, 3) 
(  II  )+HX(  11,3) 

(  II  )-tJ9.0*HX(  1  I,  2) 
(  II  )-5.C*HX(  11,2  ) 

(  II  )  +55.0+HX(  lid) 
dl  )  +  19.0*HX(  lid) 
1)+  (DT*PSUM( II)  )  /2^.0 
1), YP(2),YP(3), YP( ^),H) 

(  II  )+9.0*H{ II ) 
1)+(DT*CSUM( II) )/2^.0 
.O.DR.YCd)  .GT.l  .2)  GO  TO 
1),  YC{2),  YC(  3),YC(<.),H) 


89 
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220 


ERRQR 
6120  = 
DO  22 
DIFF- 
EL2D» 
CONTI 
£120  = 
YC3«A 
IF(El 
+NSTcP 
IFiYC 
IF(EL 
IF(EL 
IF(DT 
IF(EL 

iFcrc 

SHIFT 
DO  2  3 
DO  23 
HX(I* 
XX{  I, 
CONTI 
WRITE 
DO  24 
HX{  II 
XX{  II 
IF(YC 
N  S  T  E  P 
60  TO 
WRITE 
WRITE 
FORMA 
GO  TO 
WRITE 
FORMA 
+7X>*S 
+7X>*A 
GO  TO 
691  WRITE 
3C0  FGR".A 
+/7X, * 
♦/7X,* 
89Z  CONTI 
WRITE 
WRITE 
FORMA 
CALL 
CONTI 
CALL 
CALL 
END 


23 


24 


C 

88 
281 

S9 
290 


301 


90 


CH 
0.0 
0  I 
YC( 

EL2 
NUE 
SjP 
3S( 
2D. 
,ND 
3.G 
2D. 
20. 
.LT 
2D. 
3.G 
S  0 

J  = 

I' 
7-J 
7-J 
NUE 
(10 

II 
,1) 
>1) 
3.G 
«NS 

16 
S  V 

(D, 

T(  / 
89 
(6> 
T(  / 
-** 
FIE 
89 
(6, 
T(7 
NST 
YC3 
NUE 
(10 
(bt 
T(7 
SUV 
NUE 
GRA 
EXI 


ECK 

I-1#A 

II )-YP{II ) 

D+DIFF*bIFF*WGT( II) 


T(E 
YC  ( 
LT. 
0U3 
£.  V 
LT. 
GT. 
.DT 
GT. 
E.  V 
ATA 
1>5 
1*4 
)«H 
)»X 


DOUBLE{XX,HX,YC,H,DT/DO.ME/ 


L2D) 
3)  ) 

SEP)  CALL 

,  0  T  M  A  X  ) 

MAX)  GO  TO  391 

SEP)  GO  TO  161 

GEP)  CALL  HALVE {XX,HX,DT,DONE>NHALV) 

MIN)  GO  TO  88 

GEP)  GO  TO  161 

MAX)  GO  TO  691 


X( I>6-J) 
X( I,t-J ) 


)  YC(1)» YC{2)> YC(3)*YC(4) 

«1»4 

=  H(II  ) 

=  YC{  II) 

E.VMAX)  GO  TO  891 

TEP  +  1 

1 

ALUES  OF  SOLUTION  IF  TRAJECTORY  IS  LOST 

281)  DT,DTMIN 

7X>*ST0P!  DT«*/E20.10,2X,*DTMIN«*,E20.10) 

2 

290)L>YC(1)»YC(2)*YC(3),NSTEP 

7X, ♦ORBIT  T^n^+LOST  AT*/ 

E20.10f2X,*U«*»E20.10*2X>*V«**E20.10/ 

R  *>17>*STEPS*/  ) 

2 

300)  L,NVMAX( L) *NSTEP^YC3 

X,*G-^bIT  ♦*  I3>2X>*INTEGkATED  *  >  17 ,  *    TIMES  AROUND  TORUS*» 

EP«**  11^* 

««>E20.10) 

)  DGNE^DONE/ DONE/DONE 

301)  NHALV*NDOUB 
X,*NHALV«*,I7,5X>*ND0UB«*/I7) 
RZ(D0NE,L  ) 

FIX(DONE) 
T 
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C 

c 

c 
c 
c 
c 
c 
c 


SUB 
THI 

URE 
XMA 


QCH 

XLE 

XBO 

CCM 

WR.E 

XMA 

OCH 

XBO 

XL« 

IF{ 

XLE 

RLA 

RET 

CON 

AXM 

Rl 

RLA 

IF( 

RET 

END 


RDUTI 
S  SUB 
AL  IS 
SS  IS 
IF  XM 
IF  XI 
APGE 
NGTH 
IS  I 

AL«0. 
SS«1. 
ARGE" 
=  1.0 
XLENG 
XL.GT 
MGTH( 
R  ,1 »  0  . 
URN 
TINUE 
»AB 
-SQ 
RP  =  R1 
X.'-.ASS 
URN 


COMLARM(RLARM*L) 

TINE  COhPOTES  THE  GYRO-R AD lUS* RL ARM  . 

PUT  IN  MEGA  ELECT^-ON-VOLTS. 

PUT  IN  ATCMIC  MASS  UNITS. 

GT  CO, PARTICLE  IS  AN  ION. 

LT  O.O^PARTICLE  IS  AN  ELECTRON. 
INPUT  IN  UNITS  OF  AN  ELECTRON  CHARGE 
INPUT  IN  CENTIMETERS. 
T  IN  KILO-GAUSS. 
71/XLENGTh(5) 


NE 

ROU 
IN 
IN 

ASS 

ASS 

IS 

IS 

NPU 

ARl 

05 

0 

1.0 


TH(L) 
.0.0)  GO 
L)-1.E8 
0 


TO  1 


S(XMASS) 

RT(2.*AXM*WREAL)/QCHAPGE/XB0/XLENGTH(L) 
♦1.0218528E2 
.LT.0.0)RLAPM«R1*2.36^^&38 


SUBROUTINE  HALVE ( X X> H X, OT, DON E , NH AL V ) 

THIS  SUBROUTINE  HALVES  THE  INTEGRATION  STEP  SIZE. 

DIMENSION  KXl'f,t),HXlii,b),Hl^) 

NHALV=NHAL V+1 

DO  1  11=1,^ 

Y1"0.25*XX(II,2)+0.75*XX(II»1)-0,25*DT*HX(II,1) 

Y2-0.25*XX(II,3)+0.7S*XX(II*2)-0.25*DT«HX(II>2) 

XX{ II#6)=D0NE 

XX( II,5)«XX(II,3) 

XX( II>3)=XX(I 1,2) 

XX( 11,2) «Y1 

XX( II,^)-Y2 

HX( II,5)-HX(I  1,3)  ' 


L 
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HXdlf  3)«HX(II,2) 

CALL  EVAL(XX(1»2)>XX(2»2)»XX(3>2)>XX(^*2),H) 

DO  2  II-l*^ 

HX(II*2)«H(  II  ) 

CALL  EVAL(XX(1,^),XX(2*^)>XX{3>^)>XX(4/^),H) 

DO  3  II«1*A 

HX( II,^)-h(  II  ) 

DT»0. 5*DT 

RETURN 

END 


10 


11 


12 


SUBROUTINE  DOUBLE (XX/HX>YC>H*DT^ DON E*NSTEP> 

NDOUB,DTMAX) 
THIS  SUBROUTINE  DOUBLES  THE  INTEGRATION  STEP 
DIMENSION  XX( -^^6)  >HX  (  ^,  0)  ,  YC( -^j  >H(^) 
DIFF«A3S(XX(l,6)-DaNE) 
IF(DIFF.LT.l.O)  GO  TG  10 
IF{DT.GE.DTNAX)  GO  TO  10 
DO  1  ll^lf'f 
XX(  IIU)=YC(  II) 
XX(  II>3)-XX(I  I*'.) 
XX(  11/^  )«XX(I  I>6) 
XX(  11*6) =DONt 
XX(  Ilf 5)  -DONE 
HX(  II,1).H(  II  ) 
HX(II>3)  ='HX(  11/^) 
HX{  II>^)«HX(I  1*6) 
HX(  II*5)-DQNE 
HX{  II>6) «DONE 

WRITE (10)  YC(1)*YC(2)*YC(3)>YC(A) 
NSTEP»NSTEP+1 
DT    »OT+DT 
hDOUB»NDDUB+l 
RETURN 
CONTINUE 
DO  11  J-1/5 
DO  11  Ilxl>4 
HX{  II, 7-J) =HX(  II*6-J  ) 
XX(  II, 7-J  )  «XX(  II,o-J  ) 
WRITECIO)  YC(1)*YC(2),YC(3),YC(^) 
DO  12  11=1,^ 
HX(  11,1 )»H(  II) 
XX(II*1) -YC (I  I) 
NSTEP«NSTEP+1 
RETURN 
END 


SIZE 
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C 
C 


c 
c 


SUB 
THI 
(S* 
COM 
COM 

COM 

COM 

COM 

UPl 

VPl 

UP2 

VP2 

UP3 

VP3 

CU 

SU 

CV 

SV 

CU2 

SU2 

CU3 

SU3 

CU3 

SU3 

CUV 

SUV 

CU2 

SU2 

CU3 

SU3 

CSU 

SIN 

CSV 

SIN 

KRA 

KZA 

CAL 

CAL 

CAL 

RAD 

AN 

PN 

RAD 

RAD 

ROB 


(R/Z)  COORDINATES  OF  A  POINT 


ROUTINE  COTRAN(S»U»V>P/Z) 

S  SUBROUTINE  COMPUTES  THE 

U>  V). 

M0N/VARl/?I2*EPiGLZ*ALF 

MON/VARIA/RCC (6,10>10)>RCS(6,10>10),RSC(6*10,10)*RSS{6*10»10)/ 

IROfJROiKRO^lDR/JDRiKDR 
MDN/VAR15/CSU2iSIf;U*CSV2,SINV 

MCN/VAR19/RAC(26)>RAS{26)f2AC(26)fZAS(26)/KRAX,KZAX 
M0N/VAK20/DEL0>DEL1*0EL2*0EL3/DEL10>DEL20,0EL30*DEL22,DEL33 

«PI2*U 

-PI2*V 

«2.0*UP1 

«2.0*VP1 

»3.0*UP1 

«3.0*VP1 

»CQS(UP1) 

«SIN(UP1) 

»C0S(VP1) 

«  S  I N  (  V  P  1  ) 

»CU*CU-SU*SU 

»2.*CU*SU 

»CU2*CJ-SU2*SU 

»SU2*CU+CU2*SU 

V  "CUB+CV+SUB^SV 

V  =SU3*CV-CU3*SV 
"CU*CV+SU*SV 
•SU*CV-CU*SV 
=CUV«CUV-SUV*5UV 
=2.*CUV*SUV 
«CU2V2*CUV-SU2V2*SUV 
«SU2V2*CUV+CU2V2*SUV 
=2.0*CU 
-SU 
»2. 0*CV 

=>sv 

X1»KRAX-1 

X1«KZAX-1 

L  F0URIER(RAC*RAS>KRAX1>RAX>CSV2»SIKV) 

L  F0URIER(Z&C>ZAS>K2AX1,ZAX,CSV2>SINV) 

L  FSUB (RCC»RCS»RSC>RSS*IRO*JRO»KR0,S/RD>I0RfJDR^KDR) 

=l.-0EL0*CV+DELi0*CU+DEL20*CU2+DEL3C*CU3-0EL3*CU3V 
+DEL22*CU2V2+DEL33*CU3V3 

•3.0 

»-1.0/AN 
1   -1.0+AN*DEL3*CU3V 

■RAD1**PN 

«RAD*CU+DEL1*CV-DEL2*CUV 


V2 

V2 

V3 

V3 

2 

U 

2 

V 
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ZOB    »RAD*SU+DEL1*SV+0EL2*SUV 

R      ■RAX+S*RO*(ROB-RAX) 

Z      «ZAX*S*RO*(ZaB-ZAX) 

RETURN 

END 


SUBROUTINE  F0URIER(X>Y>K0,A,CSW2/SIiMW) 
C   THIS  SUBROUTINE  COMPUTES  A  GNE-D IMEN S ICN AL  FOURIER  SERIES 
C   BY  A  METHOD   OF  J. M. WATT, 

C    'A  NOTE  ON  THE  EVALUATION  OF  TRIGONOMETRIC  SERIES'* 
C     COMPUTER  JOURNAL, VOL. 1>»^  (1959). 

DIMENSION  X(1),Y(1) 

TRX1«0.0 

TRY1»0.0 

TRX«X(K0+1) 

TRY»Y(K0+1) 

DO  10  K^'l^KO 

TX«TRX*CSW2-TRXl+X(K0+i-K) 

TY«TRY*CSW2-TRY1+Y( KC+l-K) 

TRX1»TRX 

TRY1«TRY 

TRX  =TX 

TRY  "TY 
10  CONTINUE 

AX»TRX-TRX1*0. 5*CSW2 

A  »TRY1*SINW 

A  -A+AX 

RETURN 

END 


SUBROUTINE  PO L Y ( A, I  I  * S» VAL ) 
C   THIS  SUBROUTINE  COMPUTES  THE  VALUE 
C   DEGREE  II  BY  THE  HJRNER  METHOD. 

DIMENSION  A(l) 

VAL  "=A(  II  +  l) 

IF(II.EO.C)  -RETURN 

DO  10  I«l,  II 

VAL  »VAL'S+A(  II  +  l-I  ) 


OF  A  POLYNOMIAL  OF 


-SA- 


ID 


CONTINUE 

RETURN 

END 


C 
C 

c 
c 
c 


601 


602 

3 

5 


SU 

TH 
OR 
OR 
NU 
NY 
CO 
CA 
DO 
WR 
FO 
DO 
CA 
IF 
WR 
FC 
CO 
CA 
CO 
CO 
CA 
RE 
EN 


BROUTI 
IS  SUB 
BITS. 
BPLT3. 
MBER  I 
U  CALC 
MIDN/V 
LL  PLO 

3  LL- 
ITE(6f 
RMAT(  / 

<»  L  » 
LL  ORB 
( XPER. 
ITE  ib, 
RMAK  / 
NTINUE 
LL  FRA 
NTINUE 
NTINUE 
LL  PLO 
TURN 
D 


NE  &RAFIX(D0NE) 

ROUTINE  SUPERVISES  THE  PLOTTING  OF  THE  LINES/ 
THE  ACTUAL  PLOTTING  IS  DONE  BY  SUBROUTINE 
THE  EXTERNALS  PL OT SB L> F RAME ^ P L GT> SYMBOL , 
N  THIS  SUBROUTINE  AND  IN  OR3PLT3  REFER  TO  THE 
OMP  PLOTTER  P AC KAGE* S COOP  . 
Ak21/ VPL { ^)f NORfl/XPEP 
TS3L(60/21H  ORBIT  PLOTS  ) 

601)  VPL(LL) 

7X,*SECTI0N  AT  V«*>F6.3) 
1/NQRB 

PLT3(D0hE,L*VPL(LL)*R0T,XPER) 
GT.2.  )  GO  TO  't 

602)  L*ROT 

lOX, ♦ROTATION  NUf^.BER  OF  LINE  *>  I  3*  *  «♦>  E  20  .  10  ) 


T(0. 0*0. 0/999) 


SUBROUTINE  SU VRZ ( DON t * L ) 
C   THIS  SUBROUTINE  SUPERVISES  THE  COORDINATE  CHANGE  OF  VARIABLES 
C   OF  THE  SOLUTION  (SiU>V)  TO  (R,Z). 

REWIND  10 

REWIND  L 
10  KEAD(IO)  S,U>V,RHO 

IF ( AbS(S-CGNE ) .LE .1.0)  GO  TO  20 

CALL  CGTRAN( SiU> V>R» 2 ) 

WRITE  (L)  k,I,\/ 

GO  TO  10 
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20 


CONTINUE 
WRITE(L) 
RETURN 
END 


DONE*DCNE>DONE 


10 


SUBROUTINE  0RBPLT3(DDNE,L>VPLANE*R0T»X2) 

THIS  SUBROUTINE  PRODUCES  PLANAR  CROSS  SECTIONAL 

PLOTS  OF  THE  LINE/CRBir  IN  THE  PLANE  V«VPLANE. 

C0MMCN/VARl/PI2,EP,aLZ>ALF 

DIMENSION  A0RI&(2)>R(2)>Z(2)>V(2) 

REWIND  L 

READ( L)  R (1)>  2(1)> V(l) 

IF(AES(R(1 )-DONt ) .LT. 0.0001)  GO  TO  30 

AORIG(l) »5.5 

A0RIG(2)  >5 

SCAL 

NPER 

IND 

XIND 


5 

•3.8 
«1 
«1 
«A3S( X2-QLZ  ) 


IF{ XIND. LE,0. 00001 )IND=2 

SUMTH    »U.0 

CALL  PLOTC ADP  IG(1 )> AQkIG(2 )*-3) 

IF(L.GT.l)  GO  TO  10 

SMA«0.0 

UMA=0. 0 

VMA«VPLANE 

CALL  COTRANCSrtA^UMA*  VfiA/RMA»  ZMA  ) 

RMA»RMA*SCAL 

ZMA«ZKA*SCAL 

CALL  PLOT(  O.Of-4.5^3) 

PLOT(  0.0>  4.^,2) 

PLOTC-'t.S/  0.0»3) 

PLaT(  ^.^,     0.0/2) 

PL3T(  t^.b,     0, 


CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CONTI 

FEAD( 


0/3) 

PL0T(RMA,ZHA,3) 

SYMeOL(PMA/Zi'>lA/0.07,lHX/0.0/l) 

PL0T(PMA>ZMA,3) 

NUE 

L)  ^KZ),1K1)  fMKZ) 
VA3S-A9S ( V( 1)-V(2 ) ) 

1F( AcS( R(2 )-DGNE ) .LT.O. 0001)     GO    TO    30 
IF( VABS.GE.1.0)     GO    TO    20 
VMAX«AfiAXl(  V(  1)*  V(2)  ) 
VMAX«AMOD( VMAX,X2) 
VMAX»X2+VMAX 


-se- 


ll 


12 


13 


131 


31 


20 
600 

30 


VMAX 

VABS 

VABS 

VABS 

VMIN 

IF(V 

GO  T 

CUOT 

RAPP 

ZAPP 

IF(V 

GO  T 

RAPP 

ZAPP 

CONT 

RAPP 

ZAPP 

OIFR 

DIFZ 

IF(N 

TTHl 

TTH2 

THET 

THET 

Sur,T 

NPER 

OIFR 

DIFZ 

CALL 

CALL 

CALL 

CONT 

Rd) 

Z(l) 

V(l) 

GO  T 

kRIT 

FORM 

RETU 

CONT 

ROT 

IF(  I 

VALU 

CALL 

CALL 

CALL 

CALL 

RETU 

END 


-AMOD 
■AMOD 
«X2  +  V 
»  AMOD 
■  VMAX 
MAX.G 
0  31 
»(  VPL 
•R(l  ) 

«:Z(1) 

(1)  .6 

0  13 

»R(2) 

-Z(2) 

INUE 
=  PAPP 
«ZAPP 
2«RAP 
2  =  ZAP 
PER.E 
«OIF 
=  DIF 
A«TTH 
A«ATA 
H»SUM 
«=  N  P  E  R 
1«DIF 
1»DIF 
PLOT 
SYMB 
PLOT 
INUE 
=  R(2) 
'HZ) 
«V(2) 
0  10 
E  (6^6 
AT(  /2 
RN 
INUE 

-SUM 

ND.EQ 

E-L+O 

PLOT 

NUMB 

PLOT 

PLOT 

RN 


( VMAX,X2) 

(  VA8S*X2 ) 

ABS 

( VA8S»X2) 

-VABS 

E. V? LANE. AND. VM IN. LE.VP LANE  )  GO 

ANE-VMIN ) /( VMAX-VMIN) 
+  &UOT*(R (2)-R(l)  ) 
+  QU0T«(Z(2)-Z(1)  ) 
T,V( 2)  )  GO  TO  12 

+  QU0T*(R(1)-R(2)  ) 
+  QUOT»(Z (1 )-Z(2)  ) 


TO  11 


*S 
»S 

P- 
P- 
Q. 
Z2 
R2 
1/ 
N( 
TH 
+  1 
R2 
Z2 
(R 
OL 
(R 


CAL 
CAL 

PMA 

ZMA 

1) 

♦DI 

*DI 

TTH 

THE 

+  Th 


GO  TO  131 

FR1-DIFZ1*DIFR2 

FR1+DIFZ1*DIFZ2 

2 

TA) 

ETA 


APP>ZAPPj3) 

(RAPP,ZAPP>0.07*1H.,0.C>1) 
APP>ZAPPf  3) 


00) 

X>*IN    CR3PLT3/THE    STEP     IN    V    IS    TOO    BIG*) 


TH/FLQAT (NPER  ) 

.1)  R0T=R0T/PI2*GLZ 

.001 

(RAPPiZAPP^  3) 

ER(RAPP, ZAPPED. 1<»,VALUE>0.0/-1) 

(RAPP>ZAPP,3) 

( -A3 RIG(1)/-A0R 16(2)^-3) 
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SUBRO 
THIS 

COMHG 
COMMQ 

•f 

+ 
C0M.10 

♦  SSI 
+   CC2 

♦  CC3 
CDMMD 
COMMC 

COMMO 

M«83 
R  E  W  I N 
READC 
READ( 

+  ALF 
00  1 

.  READ( 
KZAX« 
DO  2 

:  REAO( 
REAO( 
CALL 
R  E  A  D  ( 
CALL 
R£AD( 
CALL 
PEAD( 
CALL 
READ( 
CALL 
READ( 
CALL 
READC 
CALL 
RETUR 
END 


UTINE  RE 
SUBRCUTI 
N/VARl/P 
N/VARIO/ 


N/VARll/ 

K(6>10>1 
K{6, 10,1 
Kib, 10,1 
N/VAR13/ 
N/VAR131 
N/ VARl^/ 

N/VAR19/ 
N/VAR20/ 

D  rt 

M)  iRn,i 

."JDELCD 


ADl 

NE  S 

12>t 

CCK 

CC2( 

CC3( 

CCIK 

C), 

0),C 

0),C 

L"-.  8  , 

/irK 

RCC  ( 

IRQ, 
KAC( 
DELO 


UPERVISE 
P, 3L2,AL 
6,10, 10) 
6, 10, 10) 
6,10,10) 
(6, 10, 10 

S2K(6, 10 
S3K(6,  lu 

J  M  E  ,  K  -i  B  , 

6,10, 10) 
JRD,  KRO, 
26),RAS( 
,DEL1, OE 


Mi,  IT,  IT,  IT, 

EL1,DEL2,0EL 


S  THE  INPUT  OF  DATA  FROM  TAPE83. 

F 

,CS1( 6, 10, 10), 501(6,10,10), 551(6,10,10), 

,C52(6,10,10),5C2(6,10,10),5S2(6,10,iO), 

,CS3(6,10,10),SC3(6,10,10),SS3(6,iO,10) 

),CSlK(6,10,iC),5ClK(6,10,10), 

,10),SC2K(6,10,10),SS2K(6,10,10), 

, 10),SC3K(6,1C,10),SS3K(6,10,10) 

IDB,JDe,KD3 

, IDK, JDK,kDK 

,RCS(6,10,10),P.SC(6,10,10),RSS(6,10,10), 

IDR, JDP,KDR 

26),ZAC(26),ZAS(26),KPAX,KZAX 

L2,DEL3,DEL10,DEL20,DEL3C,DEL22,OEL33 


IT,  IT, t P,QLZ,KRAX 
3,DEL10,DEL20,DEL30,DtL22,DEL33,NVAC, 


K-1, 
M  )  R 
KRAX 
K=l, 
M)  Z 
M)  I 
CREA 
M)  I 
CREA 
«)  I 
CREA 
«)  I 
CREA 
M)  I 
CREA 
M)  I 
CREA 
f-.)  I 
CREA 
N 


KRAX 

AC  (K  ),RAS (K  ) 


f<ZAX 
AC  (K 
RO,  J 
D(M, 
M3,  J 
0(M, 
KB,  J 
D(r,, 
M3  1, 
L(r., 
KK,  J 
D(H, 
MK,  J 
D(M, 

hKi, 

D{M, 


),ZA 

RO,K 
RCC, 
MB,K 
CC2, 
f,3,K 
CC3, 
J  :-^  c  , 
CCl, 

riK,K 

CC2K 

Mk,  k 
CC3K 
Jf-.K, 
CCIK 


S{K) 

RO 

RCS,RSC, 

MB 

C52,SC2, 

MB 

CS3,SC3, 

KM  3 

CSl, SCI, 

MK 

,CS2K, SC 

MK 

,CS3K,SC 

KMK 

,CS1K, SC 


RSS,IRU,JkO,KPO,IDR,JDR,KDR) 

552,  1Mb, JM3,KMB, lOB, JDB,KDB) 

553,  1MB,  JMB,KMB, IDB, JDB,KDB) 
SSl,IMEl,JMB,KMB,IDB,JDb,KD3) 
2K,5S2K,IMK,JMK,KMK,IDK,JDK,KDK) 
3K,5S3K,IMK,JMK,KMK,IDK,JDK,KDK) 
lK,SSlK,IMKl,JM<,KMt<,IOK,JDK,KDK) 
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SUBROUTINE    CREAD(M>CC*CS*SC>SS/IUJJ»KK,I1,J1>K1) 

DIMENSION     CC(Il,Jl/Kl),CS(Il»Jl>Kl)>SC(Il>Jl^Kl),SS(ll*Ji>Kl) 

CALL     READ2(M>CC,II>JJiKK>Il,Jl/Kl) 

CALL     READ2(M,CS^II/JJ*KK,I1,J1,K1) 

CALL    READ2(M,SC*II/JJ/kK>I1*J1>K1) 

CALL  READ2(M,SS/II*JJf<K,Il,Jl,Kl) 

RETURN 

END 


SUBROUTINE  READ2(M,F*IL*JLfKL>ID*JD>KD) 

THIS  SUBROUTINE  ACTUALLY  READS  TAPEMWhERE  M^BS 

DIMENSION  F(ID*JD/KO) 

DO  1  I»1>IL 

DO  1  J»li  JL 

DO  1  K«1*KL 

REAO(M)  F{I,J,K) 

RETURN 

END 


SUBROUTINE  FSUBCCC>CS>SC>SS,II/JJfKK,S,T»Ii>Jl,Kl) 
C   THIS  SUBROUTINE  CONPUTES  A  SERIES  CONSISTING  OF 
C   POLYNOMIALS  IN  S^SINES  AND  COSINES  IN  L  AND  V. 

COMMCN/VAR15/CSU2^SINU,CSV2>SINV 

C0MMDN/VAR18/P(li>)>0(15)^PP(15)iQQ(l^)fRR(6) 

DIMENSION  CC(I1>J1^K1),CS(I1,J1>K1),SC{I1>J1^K1)>SS(I1*J1/Kl) 

IIl-II-l 

JJ1«JJ-1 

KK1«KK-1 

DO  10  I«1,II 

DO  20  J«1*JJ 

DO  30  K'l^KK 

P(K)»CC (IiJ»K) 

QIK) «CS( 1/ J><) 
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30  CDSTINUE 

CALL  FaURIER(P,Q>KKl,PP(J)*CSV2>SINV) 

OD  AC  K-1>KK 

P(K)  -SC (  I/J»K  ) 

0( K)»SS (  If J>K) 
AO  CONTINUE 

CALL  FQUPIER(P»C*KK1>QC1{J),CSV^>SINV) 
20  CONTINUE 

CALL  F0UPIER{PP>QQ/JJ1,RR(1)^CSU2>SINU) 
10  CONTINUE 

CALL  POLY{RR>  II1>S> VAL) 

T-VAL 

RETUkN 

END 


C 
C 


SU 
TH 
SY 
CD 
CO 


+ 

+ 

( 

+ 

+ 


CO 


CO 
CO 
CO 
CO 
01 
OK 
DK 
DK, 
UP 
VP 
CS 
SI 

cs 

S  I 

I,-; 

CA 
CA 
CA 
IF 


BRGUT 
IS  SU 
STEM 

M  M  G  N  / 
MM  ON/ 


MMGN/ 
SS1K( 
CC2K{ 
CC3K( 
MMCN/ 
M  M  0  N  / 
M  M  0  N  / 
MMuN/ 
MENS  I 
GS«0. 
GU>=0. 
GV  =  0. 
1  -PI 
1  -PI 
U2=2. 
NO- SI 
V2-2. 
NV«S  I 
3 1  -  I  M 
LL  FS 
LL  FS 
LL  FS 
(AbS( 


I  N  E  E  V 
BRuUTI 
OF  ODE 
V  A  R  1  /  P 
VAkIO/ 


VA^ll/ 

6>10,1 

6,  10^1 

t, 10»1 

VAr 13/ 

VAR131 

V AR13/ 

V  A  P  1  7  / 

J  N  ri  (  t 

0 

0 

0 

2*U 

2*  V 

0*CGS( 

'i  (LP!  ) 

0  *  C  0  S  ( 

\  (  V  P  1  ) 

6  +  1 

U3{CC1 

U3{CC2 

UB(CC3 

RLAKM ) 


AL(Si 

NE   C 
S 

12, EP 
CCi  (6 
C  C  2  (  6 
CC3(  b 
CCIK  ( 
0), 
0)>CS 
0)  >CS 
I  M  B  ,  J 
/  I ,''.  K  , 
CSU2, 
RL  ARM 
) 


UMPUTES  THE 


RIGHT  HAND  SIDE  OF  THE 


,CLZ, ALF 

,  iGilO)>CSl{6,10,lG),SCl(6»10,lC),SSl(b>lO>lG)> 

>lC/10)KS2(6flO,lC)>SC2(6>10>10),SS2(6,lC,10). 

>lC>10)iCS3(6,10,10),SC3(5,10,10),SS3(t,,lGfl0) 

0,10>10)/CSlK(6flO,lC)>SClK(6>10,10), 

2i<{6,lCflO)>SC2K(6»lG,lC)>SS2K(t,  lO^lD), 
JK.  (6»10/10)>SC3K(6,1C/10),SS3K(6*10,10) 

M  c  ,  K  M  B  >  I  [;  b  ,  J  D  b  f  K  D  B 
jrK^KMK, IDK/ JDKfKDK 
SINU,CSV2/S INV 


UPl  ) 
VPl  ) 


,CSl>SCl,SSl,IMBl/JMB,KMB,S*DBGSiIDB>JDB,KDB) 
>CS2,SC2»SS2,IMbfJMB,KMBfS,DBGU#lDB>JDB»KD8) 
>CS3,SC3>SS3,IMB,JMB,KMB»S>0BGV*I0B>JDB>KDB) 
.LT.l.E-6)  GO  TO  1 
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IMK1»I 
CALL  F 
CALL  F 
CALL    F 

CCNTIN 
H{1)>»D 
H(2)-D 
H(3)0 

RETURN 
END 


[MK  +  1 

-SU3(CC1K,CS1K>SC1K»SS1K,IMK1, JMK,KMK,S*DKGS, IDK, JDK,KDK) 

■SUB(CC2K,CS2K,SC2K,SS2K^IMK#JMK,KMK,S/DKGU*IDK,JDK,K0K) 

■SUB(CC3K,CS3K,SC3K*SS3K,  I  «K,  J  rlK,  KMK*  S*  OKG  V>  IDK,  JDK,KDK) 

<UE 

)3GS+RLARM*DKGS 


3GS+RLARM+DKGS 
^GU+RLAR.-^+DKGU 
6GV+RLARM*DK6V 
.0 
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LISTING  OF  THE  CODE  COEFF 


PROGRAM  C0EFF(INPUT>0UTPUT,TAPE6«0UTPUT>TAPE7,TAPE82* 
♦  TAPE83/TAPE8>TAPE9/TAPE10) 

C   THIS  IS  THE  DRIVER  ROUTINE  THAT  CALLS  THE  OTHER  SUBROUTINES 
C   TAPE?  IS  THE  INPUT  DATA  FROM  THE  EQUILIBRIUM  CODE 
C   TAPE63  IS  THE  OUTPUT  TAPE  FOR  LIN0R3  CODE 
C   TAPES  3>9, 10*82  APE  FOR  TEMPORARY  STORAGE 

IREAD«1 

CALL  READKIREAD) 

CALL  UNIV(IREAD) 

CALL  CNSN 

CALL  SUBAX 

CALL  SU3R0 

CALL  SUBPSI 

IS-1 

CALL  DBGRADS(IS) 

IREAD-3 

CALL  UNIVdREAD) 

CALL  CNSN 

CALL  SU3HJK 

CALL  EXIT 

END 


SUBRCUTIN 
THIS  SUBR 
INITIALIZ 
BY  THE  DA 
IMAXRO 

JMXRO 

KMXRO 

I M  A  X  P  S  I 

JMXPSI 

KMXPSI 

KAX 


E  UNIVdRE 
DUTINE  DEF 
ATION  OF  T 


TA 
-1 

-1 
-1 
-1 

-1 
-1 
IS 


IDEG-1 


JMAX-1 


ATEME 
DEGR 
IN  T 
NO. 
NO. 
DEGR 
DBGR 
NO.  , 
NO. 
NO, 
COOR 
DEGR 
DKGR 
NO. 


AD) 

INES  PARAM 
HE  RUN. THE 
NTS  BELOW. 
EE  OF  POLY 
HE  MAPPING 
OF  HAPMONI 
OF  HARMONI 
EE  OF  POLY 
ADU  AND  DB 
OF  HARMCNI 
GF  HAkMCNI 
OF  HARMONI 
DINATE  FUN 
EE  OF  POLY 
ADU  AND  OK 
OF  HARMONI 


ETERS  FUR  THE  PROPER 
Sb  PARAMETERS  ARE  DEFINED 
THEIR  MEANINGS  ARE  AS  FOLLOWS: 
NOMIAL  APPROXIMATION  OF  RO 

• 

CS  IN  U  OF  RO.   • 

CS  IN  V  OF  RO. 

NOMIAL  APPROXIMATION  OF 

GRADV. 

CS  IN  U  OF  OBGRADU  AND  DBGRaDV. 

CS  IN  V  OF  ObGRADU  AND  DBGRADV. 

CS  IN  V  OF  RAX  AND  ZAX,THE(R*Z) 

CTiCNS  OF  THE  MAGNETIC  AXIS. 

NOMIAL  APPROXIMATION  OF 

GRADV 

CS  IN  U  OF  DKGRADU  AND  DKGRADV. 
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C         KMAX- 
C    (SI/UI,VI 

c 

c        w 

C  VK 

C  (IR^JRiKR. 

C 

G  (IA,JA,KA 

c 

C    (IB*JB,KB 

C 

C    (ICJCfKC 

c 

C    NNI*NNJ*N 
C 

COMMON/VAR 
COMMQN/VAR 
CGM-iON/VAR 
COM'iGN/VAR 
COMMON/VAR 

+ 
COMMON/VAR 

+ 
COMMON/VAR 

+ 
COMMON/VAR 

+ 
COMMON/VAR 

+ 
CCMMGN/VAR 
COMMON/VAR 

+ 
CCMMCN/VAR 
COMMON/VAR 
COMMON/VAR 

+  KAXIS* 
DIMENSION 
DATA  IMAXR 
DATA  IMAXP 
DATA  KAXIS 
DATA  IDEG^ 
DATA  SI^UI 
DATA  W»VK, 
DATA  NNI*N 
DATA  IR/JR 
DATA  IA,JA 
DATA  IB»JB 
DATA  ICJC 
IF{  IREAD.E 
IF( IREAD.E 
100  CONTINUE 

ITEXT(1)"6 


1   ■  NO.  OF  HARMONI 
)   -  (S/U/V)  COOROI 
THE  PARTICLE. 

•  ENERGY  OF  THE 

-  INITIAL  VALUE 

)   »  DIMENSION  OF  T 

COEFFICIENTS  0 

)   =  DIMENSION  OF  T 

COEFFICIENTS  0 

)   «  DIMENSION  OF  T 

COEFFICIENTS  C 

)   »  DIMENSION  OF  T 

COEFFICIENTS  0 

NK  '    NO.  OF  S»U*V  P 

IN  THE  CQMPUTA 

1/NI,NJ,NK>NVAC*  IDE 

2/HSf HU>HV 

3/PI2>EP^QLZ* 

d/IPU>  IPV^ IPU 

11/CSU2>SINU> 

PPP( 10) 
13/IR,JR^KR,IR0>JR0 

RSC (6*10ilO)*RSS 
l-i/IB*  JB^KB>  IBU>  JBU 

SC2(6>  10*10)  /SSi: 
15/  IC> JC>KC*  IBVt  J6V 

SC3(6> 10>10)>SS3 
16/  lA/ J  A*KA* IBS* J6S 

SC1(6,1G*10)*SS1 
17/KRAX,KZAX*  RACdO 
IB/DELOi  DELI*  DEL?*D 

DEL33 
19/W*  XMU*  VK 
20/JMAX*KMAX* 
22/  IMAXRO* IMA 
JMXRQ*  KMXRO*  J 
ITEXT(5) 
0*  JMXROfKMXRO 
SI*JMXPSI*KMX 
/5/ 

JMAX*hMAX/3*b 
, VI  /O. 66657/0 
NRE  S/O.  5^-0.6 
NJ*  NNK/13*  Z^, 
/KR/6/10/10/ 
/KA/6/10,lC/ 
/KB/6/10/10/ 
,KC/6,10* 10/ 
C.l)  GO  TO  100 
Q.3)  GO  TO  300 

HCONSTANT 


CS  IN 
NATES 

PARTIC 

OF  VPA 
HE  ARR 
F  RO. 
HE  ARR 
F  OBGR 
HE  ARR 
F  D3GR 
HE  ARR 
F  DBGR 
OINTS* 
TIQN  0 
G 


V  OF  DKGR 
OF  INITIA 

LE. 

R/VPERP  0 
AYS  CONTA 

AYS  CONTA 

ADS. 

AYS  CONTA 

ADU. 

AYS  CONTA 

ADV. 

RESFECTIV 

F  DKGRADU 


ADU  AND  DKGRADV. 
L  POSITION  OF 


F  THE  PARTICLE. 

INING  THE 

INING  THE 

INING  THE 

INING  THE 

ELY, 
AND  DKGRADV. 


ALF 

1/IPVl 

CSV2*SINV*P(10),Q(10) 


*KRO>R 
(6,10/ 
,KBU*C 
(6,10> 
>KBV,C 
(6,10, 
,KBS>C 
(6/10, 
),RAS( 
EL3,DE 


CC(6,10/1 

10) 

C2(6,10,l 

10) 

03(6,10,1 

10) 

CI  (6,10,1 

10) 

lO/ZACd 

L10,DEL20 


,PP(10)/QQ(10), 

0)/RCS(6/10/10)/ 

0),CS2{6, 10,10), 

0),CS3(6,10,10), 

0),CS1(6,10,10), 

0), ZAS( 10) 
,DEL30,DbL22, 


S  0  ,  N  R  E  S 

X  R  S  I ,  N  M  A  X  , 

MXPSI,  KMXPSI 

/3,6,3/ 
PSI/3,6,3/ 

,3/ 

.0,0.0/ 

,9/ 

2-;/ 
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I 
I 
I 
I 
P 
N 
I 
N 
W 
10  F 
+ 
I 


300 


301 

302 
303 
304 

305 
306 


TEXT(2 ) 
TEXT(3) 
TEXT(<») 
TEXT(5) 

I2«2.0* 
KAX»10 
MAXVAC" 
IV»0 
RITE (6^ 
ORHAK  / 

F{NJ/2* 
PU« (NJ- 
0  TO  3 
PU«(NJ- 
F (NK/2* 
PV=(NK- 
0  TO  6 
PV» ( NK- 
ONTINUE 
PUl'IPU 
PV1»IPV 
RITE(82 
RITE(82 

ETURN 
GNTINJE 
I»KM 
J  •  fi  N  J 
K«NNK 
0»1.0/1 
I2»2.0* 
F  (  N  V  A  C  . 
S»  (1.0- 
U  -1.0/ 

V   :'1.0/ 

RES=SO+ 
I«SRES 
F  (  N  J  /  2  » 
P  U  »  (  N  J  - 
0  TD  3  3 
PU"(NJ- 
F  ( ^^  K  /  2  * 
P V- (MK- 
U  TO  30 
PV« (NK- 
DNTINUE 
PUl'IPU 
PV1»IPV 
ALL  REA 
L  «  S  I  ♦,♦  A 
DER-1 


■6HLINEAR 

-QHQ'JADRATIC 

•5HCUBIC 

-7HQUACTIC 

ACOS(-l.O) 


10)  ITEXT(IK'AXRO),ITEXT(IMAXPSI) 
/5X> A10/20HAPPROXINATIDN  FOR  RQ /  /  , 

5X,A10>29HAPPR0XIMAT lUN  FOR  B  IN  PLASMA) 
2-NJ) 1,2^1 

1)  /2 

2)  /2 

2-NK)4,5,4 

1)  /2 

2)  /2 

♦  1 
+  1 

)IMAXP0*IMAXPSI>1MAXVAC^IPU1^IFV1>NI,UIV>EP>0LZ»KAXIS 
)DELG>DELl>DtL2>DEL3^DfcL10»DEL20>DfcL30^DEL22»DEL33^ 
NVAC>  ALF 


2.0 

AC0S(-1 .0) 
GT.  20)     5  0  =  0.0 
SO)  /FLuAT(M-l) 
FLOAT (NJ ) 
FLOAT (U\ ) 
(NRES-1. )*HS 

2-NJ ) 301»302>  301 

1  )  /2 

3 

2)  /2 

2-r^K)  3C4i  30i^  30-^ 

l)/2 

6 

Z)/Z 

♦  1 
+  1 

DKIREAD) 
LF 
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CAL 

XMU 
WRI 


L  EVA 
•W/(  1 
TE(6, 


WRITE(6, 


WRI 
WRI 

600  FDR 

♦ 

■f 
+ 

601  FOR 

+ 
+ 

+ 

602  FOR 

+ 
+ 
+ 
■  ♦ 
♦ 

603  FOR 
+ 

+ 

+ 
RET 
END 


TE(6, 

TE  (6, 

MAT(  1 

5 

5 

5 

5 

5 

,1 A  T  {  / 


L(SL> 

.+VK* 

600)a 

D 

601)  I 
K 

602)  W 
603)N 
HI/// 
X,*IN 
X>*AL 
X,*DE 
X,  ♦DE 
X>*DE 
/5X,« 
5X>*I 
5Xi*I 
5X,*I 
5X,*K 

MAT( //5X,* 
5  X  ,  *  K 
5X,*I 
5X,»I 
10X>* 
5X,*I 
/5X,* 
X,*HS 

X,*IO 

X  ,  ♦  I  N 


MAT(/ 
5 

5 
5 
5 

URN 


UI*  VI*  I 

VK)  /XMO 
LZ»tP> A 
EL22*DE 
BS*JB5* 
RO*KRAX 
>  XMU> VK 
I^NJ  W<K 
5X,*0LZ 
VERSE  r 
F«=*,  FlC 
LO  '*, F 
LB  '*, F 
L30  =  *>  F 
IBS  =  »>  I 
6U=*,  15 
BV  =  **  15 
RO»**lp 
RAX=*>  I 
ENERGY 
AGNETIC 
NITIAL 
NITIAL 
S*»»£20 
NITIAL 
^a  =  ♦  *  I  5 
-**F10. 
AC=**I5 
EG=*>Ip 
ITIAL  F 


DER»XM 
03 

LF/DEL 
L33 

Kas>iB 

>.KZAX 
/SI>UI 
*H5»HU 
=**Fi5 
AJGR  R 
.5// 
10.5>2 
10.5,2 
10.5>2 
5,2X,* 
/2X*»J 
,2X,*J 
>2X>«J 
5i2X,* 
OF  PAR 
M  D  !1  E  N 
VPAR/V 
PUSITI 

.  1  a  / 1 0 

VALUE 

5>2X,* 
I 

>2X,*J 
LUX  SU 


ODB,DKS,0KU,DKV) 

0>DEL1,0£L2,DEL3*DEL10*DEL20»D£L30* 

U>JBU*KBU,IBV,JBV/KBV,IRO*JRO* 

,VI,XMODB 

>HV,NVACiID£G,JMAX,KMAX,SO 

.8// 

ADIUS  OF  TQRUS»**E20.10// 

X**DEL1  =*>F10.5,2X,*DtL2  =*,F10.5/ 

X,*DEL10=*,F10.5,2X>«DEL2  0=*,F10.5/ 

X,*DEL 22 =*>F10.5>2X>*DEL33=*iF10. 5) 

JBS=*,I5*2X*«KBS=*, 15/ 

6U  =  *>  I5*2X, ♦K6U=*, 15/ 

BV-*,  I5*2X**KBV«*f 15/ 

Ra»*>  I5,2X>*KRJ  =  «> 15/ 

KZAX=*> 15) 

TICLE=**E20.10/ 

T    '=*,E20.1C7 

PEF  P=*>E20.10/ 

ON  OF  PARTICLE:*/ 

X,*U=*>E20.1C/10X,*V»*>E20.10/ 

OF  MCDB"*>E2C. 10) 

J=** 15/2X, *NK«** 15/ 

HU-**  F10.5*2X,*HV»**F10.5/ 

MAX  =  *,  I5*2X**KMAX  =  *,  15/ 
PFACt     AT       *>F10.5) 


SUBROUTINE     READKIREAO) 
C       THIS     SUBROUTINE    READS    TAPE7    (IF     IREAD»1)>AN0    TAPE82 
C        (IF     IREAD^B).     IT     CALLS     THE     SUBROUTINE    CREAD^wHICH    CALLS 
C       THE    SUBROUTINE    READ2. 

C0M•1CN/VAR1/M,^.J,NK,NVAC>IDEG 

COhiCN/VAR3/PI2,£P,CLZ>ALF 

C0MM0N/VAR8/IPU,IPV,IPU1/IPV1 

CL.1M0N/VARi3/IR,JP,KRWKG,JRU,KRD*RCC(6>10,10)>RCS(6ilC,lU), 

♦  RSC (6>10f 10)>RSS (6>10/10) 

CDMMDN/VARl^/  IB/JB,KB/IBU>JBU*KBU*CC2(6>10>10)*CS2(6,10*10)* 

♦  SC2(6,10>10),SS2(6,10>10) 
C0MM0N/VAR15/IC,JC,KC,IBV,JBV>KBV,CC3(6,10»10)»CS3(6,10*10)* 
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SC3(6>10 

C0MM0N/VAR16/IA,JA^KA 
SC1(6^10 
C0MH0N/VAR17/KRAX*KZA 
CDMiON/VARlS/DELO^DEL 
DEL33 
.EQ.l)  GO  TD 
.EQ.3)  GO  TO 


*10, 10) 

BSfCCl(6flO»10)>C51(6/10,10), 

>10,10) 

RAS( 10),ZAC ( 10),ZAS( 10) 

3*DEL10*D£L2  0*DEL30*DEL22* 


IF (  IREAD 
IP (  IREAD 
100  CONTINUE 
REAO{ 7/1 
REA0(7/2 
READ(7»  3 
+DEL33>AL 

1  FDRMAT(5 

2  FORMAT (5 

3  FORMATd 
ALF-0.5 
WRITE(6> 

+DEL0/ DEL 

10  FDRMATCl 

+  5X*  5HNVA 

+  5X./19HNa 

+21X>3HNJ 

+  21X,  3HrilJ 

+  5  X  >  3  5  ri  I  N 

+  5X,'tHQLZ 

+F8.^/5X, 

+6H0EL20" 

+6HDEL33- 

RETURN 

300  CONTINUE 

M«82 

N«83 

REWIND  H 
REWIND  N 
REAO(M) 
WRITE (N) 
R£AO(M)D 
+  N 
WRITE (N) 
+  N 
KZAX'KR A 
DO  31  <« 
REAO( M) 
WRITE (N) 

31  CONTINUE 
DO  32  K- 
REAO(M) 
WRITE (N) 

32  CONTINUE 
READ(«) 
WRITE (N) 


)NVAC/NI,NJi  N 

)HS»HU»HV>EP> 

)DEL0>DEL1/DL 

F 

15) 

F12.  J  ) 

OFB."^  ) 

10)NVAC>NI,NJ 
1>0EL2,LEL3>D 
H1//2X,*DATA 
C  =  /  15/ 
.GF  POINTS: 
«>  I5/21X>  3HKK 
*>F12, 8/21X^3 
VERSE  MAJG?  K 
»  /  F  1  2  .  e  /  5  X , -^  H 
6H0EL2  «>Fa.^ 
,F8.^/5X,6HCE 
/F8.4 ) 


*10),SS3(6 

,  IBS> JBS/K 
>10)f SSI (6 
X,RAC(iC)> 
1»0EL2*DEL 

100 
300 


K,NTAPE 

QLZ 

L2,DEL3/DEL10,0EL20>DEL30>DEL22> 


»NK>HS>HU/ 
EL10»DEL20 
FROM    EOUIL 

NI  =  >  15/ 
=  >  15/13X^1 
HHV=*  F12 .8 
ADIUS  LiF  T 
ALF=, F12.e 
/ pXi  6HDEL3 
LSO^/FS."^/ 


HV/EP/QLZ» ALF, 
/DEL30tLEL22>DEL33 
IBRIUh    CODE:*// 


IHMESH:  HS=/F12.8/ 

/ 

OPUS  :  EP«/F12.8/ 

/5X/6H0EL0    =»F8.4/5Xf dHDELI 

"/F8.^/5X>6HDEL10=»F8.^/5X/ 
5X,6HDEL22=/ F0.^/5X> 


IRO*IBU>IT»IT>IT>IT*IT/EP/OLZfKRAX 

IRO/IBUWT>IT>IT>IT/IT/EP>QLZ/KRAX 
EL0*DEL1>DEL2>0£L3»0EL10>DEL20/DEL30>DEL22»DEL33/ 
VaC^ALF 

DELC»DELl>DEL2»DEL3>DEL10,DEL2C/DfcL3  0>DEL22/DEL33> 
VAC* ALF 
X 

1,KPAX 

RaC (K) »RAS ( K) 
RAC (K  )/RAS(K) 

liKZAX 

ZAC  (K  )  ,ZAS(K) 

ZAC  (,<),ZAS(K) 

IRO/ JR0*KR0 

IRQ*  JRD/KRO 
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CALL  CREAD(M,N,PCC>RCS*RSC>RSS*IRO/ JRO,KRO,IR, JR*KR) 

PEAD{M)  IBU/ J3U,KBU 

WRITE (N) IBU»JBU*KBU 

CALL  CREAD(M,N,CC2*CS2»SC2/SS2, IBU*JBU>KBU#IB,JB,KB) 

READ(M) ieV> JBV>KBV 

-RITE (N )  IbV,JBV>K3V 

CALL    CREAD(MWOCC3>CS3>SC3>SS3, IBV>JBV^KaV/IC*JC>KC) 

READ(r) IBS, J35#KBS 

W  R  I  T  E  ( rni  b  3  ,  J  B  S  >  K  B  S 

CALL    CREAD(H,M,CCl,CSl*SCl*SSl>IBS*JtJS,K8S>IA*JA,KA) 

RETURN 

END 


SUBROUTINE    CREAD(K,N,CC>CS,SC>SS,II*JJ,KK,I1»J1,K1) 

DIMENSION    CC(Il,JliKl),CS(Il>Jl*Kl),SC(Il,Jl,Kl)>SS(Il,Jl,t<l) 

CALL     kEAD2(M,h;,CC>II,JJ*KK,Il,Jl,Kl) 

CALL     READ2(M,N>CS>I1,JJ,KK,I1,J1,K1) 

CALL     READ2{M,N,SC>II,JJ,KK,Ii,Jl,Kl) 

CALL    READ2(H>N/SS,II/JJ>KK>I1>J1>K1) 

RETURN 

END 


SUBROUTINE  READ2(M,NjFiIL*JL,KL>ID>JD>KD) 

THIS  SUBROUTINE  READS  TAPE.^.  AND  WRITES  ONTO  TAPEN. 

DIMENSION  F(ID,JD,KD) 

DO  1  1  =  1, IL 

DO  1  J«lf  JL 

DO  1  K«l , KL 

READ(M)  F(I,J,K) 

WRITE(N)F(I,J,K) 

CONTINUE 

RETURN 

END 
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SUBROUTINE    CNSN 
C       THIS     SU3'?3UTINE    DEFINES    THE    ARRAYS    WHICH    WILL 
C       THE    FAST    FOURIER     TRANSFORM    SOBROUTINE     FFDRM. 

COMMCN/VAPl/NIiNJ,NK,NVAC»IOEG 

COMrtON/VAR3/PI2^EP^OLZ>ALF 

C0MM[]N/VAP6/CNU(50)>SnU(50)>CNV(50)>SNV(50) 

DT»PI2/FL0AT( NJ) 

ANG«C.O 

DO  1  J»l»  NJ 

CNU(  J  )«COS(ANG) 

SNU( J )«SIN(ANG) 

ANG» ANG+OT 

1  CONTINUE 
DT»PI2/FL0AT(NK) 
AN6«0.0 

DO  2  K=1,NK 
CNV(K )-COS (ANG) 
SNV(K)=SIN(ANG) 
ANG«ANG+DT 

2  CONTINUE 
RETURN 
END 


BE  USED  BY 


SUBROUTINE  SUBAX 
C   THIS  SUBROUTINE  FOURIER  ANALYZES  THE  DATA  FROM  THE 
C   THE  ECUILIBRIUM  CODE  FOR  THE  COORDINATES  (R>Z)  OF  THE 
C   MAGNETIC  AXIS. 

COMPLEX  CMAT>CALX,X 

COM-iCN/VARl/MW-JWiK^NVAC^  IDEG 

C0M^1LN/VAR6/CNU(50),SNJ(:;0)KNV(5C),SNV(30) 

C0MMDN/VAR7/C.^AT(51)^CAUX(51/51)^X(51) 

C0MM0N/VAK8/IPU^IPV^IPUl/IPVi 

C0MMDN/VAR21/Fl(52)>F2(!52)>FFl(52>52)/FF2(t.2>52)/FF3(52^52), 
+  FFFi(13> 2o>26)^FFF2(13i26>2fe) 

C0MH0N/VAR22/IMAXR0^IMAXPS WNMAX, 
+  KAXIS>JMXPO>KHXRO/JMXPSI>KMXPSI 

DIMENSION     ITEXT(2) 

PI2»2.0«AC'JS(-1.0) 

ITEXT( 1 ) "3HRAC 

ITEXT  (2  )  ^'SHRAS 
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REA0(7*1)  T1#(F1(K),K-1,NK),T1,T1>(F2(K)*K«1,NK),T1 

1  F0RMAT(^E16.8) 
WRITE(6,10) 

10  FORMAK 1H1//2X,37HF0URIER  COEFFICIENTS  OF  MAGNETIC  AXIS) 

DO  2  1-1*2 

DO  3  K«1*NK 
3  C«AT(K)»CMPLX(Fl{K)/0.0) 

CALL  FFQRM(NK,CKAT,X,CNV>SNV) 

WRITE(6*A)  ITEXT(l),ITtXT(2) 
A  F0RMAT(//5X,1HK,10X,A3*20X,A3) 

DO  5  K=1,NK 

5  ChAT(K)»CMAT(K)/FLOAT(NK) 
CMAT(NK*1)»C1PLX (0.0*0.0) 
DO  6  K'lf  IPVl 

•A«REAL(CMAT(K)  +  CKAT{Nr<+2-K)  ) 
B»-AIMAG(CMAT ( NK+ 2-K ) -C M AT ( K ) ) 
WRITE(6»7)  K-1*A*B 

7  FQRfiAK  IX,  15*  2E20.10) 
IF(K.GT.KAXIS)  GO  TO  6 
WRITE(82)  A*a 

6  CONTINUE 
IF(I.E0.2)RETURN 
DO  8  K=1,NK 

8  FKK)  «F2(K) 
ITEXTd  )*3HZAC 
ITEXT(2)=3HZAS 

2  CONTINUE 
RETUPf^ 
END 


V  THE  DATA 
FUNCTION  RD 
SUPERVISES 


SUBROUTINE  SU3P0 
C   THIS  SUBROUTINE  FOURIER  ANALYZES  IN  U  AND 
C   FRO"^  THE  EQUILIBRIUM  C3DE  FOR  THE  MAPPING 
C   IT  THEN  CALLS  THE  SUERQUTINE  SELECT  WHICH 
C   THE  APPROXIf^ATICN  IN  S. 

C0MMCN/VAR1/M*NJ*NK*NVAC*  IDEG 

C0MM0N/VAR21/AI(52)*A2(32)*AAl(52*52)*AA2{52*52)*AA3(b2*52)* 
♦  Fl( 13*26*26)*F2( 13*26*2t) 

IFS«1 

LFS  =  M 

N2»NJ+2 

N5-NK+2 

READ(7*1)(({F1(I*J*K)*K«1*N5)*J=1*N2)*I-1*NI) 
1  FORMATCtElb.e) 

IU-0 
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21 


22 
2 


IV-0 

DO  2  I-IFS>LFS 

DO  21  J«1*NJ 

00  21  K-1,NK 

AAK J»K )-Fl( I f J+1,K+1) 

CALL  CF3RM(  lU^ IV»  AAl ) 

DO  12    J=»1>N2 

DO  22  K-1,N5 

Fl(  I,  J/K)«AAl(J,t<  ) 

CONTINUE 

lND-1 

CALL  SELECK  IND,  IFS,  LFS>F1) 

RETURN 

END 


C 

c 
c 
c 


SUBR 
THIS 

POIN 
FOUR 
THE 

C  0  ."i  "I 
C  0  M  M 
CO  MM 

IFS- 
LFS" 

Nl-N 
^^  2  »  N 

N^«N 
N5«N 
READ 
FORM 
HS-1 
HU«1 
HV-1 
DO  2 
S»FL 
SALE 
DO  2 
DO  2 
F2(  I 
FKI 
CONT 
IU«0 
IV-1 


DUTI 

SUB 

TS  F 

lER 

SU3R 

CN/V 

CN/V 

CN/V 

2 

M 
J+1 
J  +  2 
K  +  1 
K+2 
{  1»  1 
AT(^ 
.0/F 
.0/F 
.0/F 
1  =  1 
CAT  ( 
=  $»* 
J-2 
r<«2 
*  J»  K 
»  J»K 
INUE 


NE  SUBPSI 

ROUTINE  COMPUTES  DBGRADU  AND  DBGRAOV  AT  THE  GRID 
ROM  THE  DATA  OF  THE  E  OU I L I  bR  I  Ur.  CODE.  THEN  IT 
ANALYZES  Ifi  U  AND  V  AND  CALLS  SELECT  AS  IN 
OUTINE  SUBRO  ABOVE. 
Af.  l/NI,NJ,NK>NVACt  lOEG 
KV  3/PI2/EP>CLZ> ALF 

AR21/A1(52)^A2(  52  )iAAl(  52/52), AA2(:>2,i2)>AA3(5202)i 
Fl(13,26,26) ,F2( 13/26,20) 


)  (  (  ( 

E16. 

LDAT 

LOAT 

LDAT 

FS,L 

I-l  ) 

ALF 

,Ni 

,N^ 

)«(  F 

)«(F 


1, N5) , J=1>N2) ,I=1,NI) 


Fl(  I,  J,K),K^ 

8  ) 

(  N  I  - 1  ) 

(NJ  ) 

(  NK) 
FS 

«H5 


1(I,J+1,K)-F1(I,J/K))/HU*SALF*2.0 
1(I,J>K)-F1(I>J,K+1))/HV*SALF*2.0 


-100- 


31 


32 
3 


71 


72 

7 


DO    3     I«IFSf LFS 

DO    31    J'lfNj 

DO    31    i<»l>NK 

AAK  J*K)xFl( I>  J+l^K  +  1) 

CALL  CFORM( lU/IV* AAl) 

DO  3  2  J"1*N2 

DO  3  2  K-1,N5 

Fll  I^  J/K)*AA1  {  J^K) 

CONTINUE 

IND»2 

CALL  SELECT(INC,IFS^LFS>F1) 

IU«1 

IV«0 

DO  7  I-IFS,LFS 

DO  71  J=1>NJ 

DO  71  K=1,NK 

AA1(J>K)«F2(I^J+1>K+1) 

CALL     CFOPM(IUWV^AAl) 

DO  72  J=«l/N2 

DO  72  K=l,  N5  ;  :■■... 

F2(I*  J^K  )«AA1  (  J*K  )     ;.  , . 

CONTINUE   . 

IND»3 

CALL  SELECT( INQ>  IFS> LFS>F2) 

RETURN 

END  '  - 


C 

c 

C 

c 
c 
c 
c 
c 
c 


SUBRCUTI 
This  SUB 

1.  DAT 
IN 
SUB 
SOU 
SIN 

2.  THE 
FOR 
TAP 

C  C  M ''  D  f>i  /  V 
C  0  M  '^  D  N  /  V 
C  0  M  M  C  N  /  V 
CuMMCN/V 

+  KAXI 
DIMENSIO 

♦JIND( 31) 
DIMENSIO 


NE  SE 
ROUTI 
A  PAS 
THE  P 
RGUTI 
AkES 
TRO; 
FINA 

E?2. 
A  K  1  /  N 
A  R  3  /  ? 
AR3/I 
AR^2/ 
S  /  J  .1 X 
N  F(l 
»KIND 
N  ITf 


LECTdNDi  IFS»LFS>F  ) 

NE  HAS  TWO  FUNCTIOriS: 

SEC  TO  IT  FRUi^.  SUBRO  AND  Sbt3PiI  ARE  ARRANGED 

ROPER  t^AY  SU  THAT  THEY  .lAY  3 1  PASSED  TO  THE 

NFS  SINTRU  AND  SIKT  FDR  PCLYNOMIAL  LEAST 

FITTING  IN  S.  DATA  FROM  SUBRC  ARE  PASSED  TO 

DATA  FRG^l  SU-JPSI  ARE  PASSED  TO  SINT. 

L  COEFFICIENTS  IN  THE  SERIES  REPRESENTATIONS 

DBGkADU  and  D8GRA0V  ARE  THEN  WRTTTEN  ONTO 

I^  NJ>  NK, NV AC*  NNDEG 

12*EP,OL2, ALF 

PUi  IPV*  IPL;1>  IFVl 

I  '•■  A  X  R  Ll *  1  .'i  A  X  P  S  I  *  N  ,••1 A  X  > 

rj*khxro>jmxpsi>kmxpsi       '    ' 

3*26,26)>FVAL (13  )* PC( iO)*Y(31), 

(  31) 

XTdC) 
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10 

11 


IND 

INO 

IND 

IMA 

NP- 

DO 

Y(  I 

JIN 

KIN 

CQN 

ITE 

ITE 

ITE 

ITE 

ITE 

NMA 

Il» 

12- 

DD 

00 

SUM 

DO 

SUfi 

SUM 

DO 

IF( 

CON 

GO 

DO 

Y(N 

JIN 

KIN 

CON 

Y(N 

JIN 

KIN 

CON 

WRI 

FOR 

DO 

J.J 

K»K 

Kl. 

IF  ( 

Jl« 

WRI 

GO 

Jl« 

WRI 

FOR 

CON 

WRI 


-1  I 
-2  I 
«3  I 
X-0 
LFS- 

1  I» 
)«0. 
D(  I  ) 
0(1) 
TINU 
XT(1 
XT(2 
XT(3 
XT(^ 
XT(5 
XI  »N 
IPUl 
IPV  + 

2  J» 

2  K  = 
«C.O 

3  I- 
=  SUM 
»SQR 
tf    h' 

5  L  -I . 
TINU 
TO  2 

6  M» 
MAXl 
D(NM 
0  (  N  « 
T  INU 
)  «SU 
D(  N) 
D(N) 
TINU 
TE(6 
•i  A  T  ( 
8  N« 
I N  D  ( 
I .'.  D  ( 
K-IP 
J.Lc 
J-IP 
TE(6 
TO  1 
J-1 
TE  (6 
MAT  ( 
TINU 
TE  (6 


S    RO    INDICATOR 

S     DBGRADU     (PLASMA)     INDICATOR 

S     DBGRAD\/     (PLASMA)     INDICATOR 

IFS  +  1 
1*31 
0 
-0 

»o 

E 

)-2HR3 

)«7HD3GRADU 

)  «7HD3GRADV 

)  "^SHSIN 

)-3HC0S 

MAX  +  1 

+  IPU1 

IPV+1 

1>I1 

1*  12 

IFS^LFS 

+  ( F{  I>  J»K) *F(  I> J>K)  ) 

T( SUM/FLCAT(NP  )  ) 

1  *  N  M  A  X 

GT.Y(N) )  GO  TO  b 


*NMAX 

N-M)  =  r  {NMAX  +  N-:'i  ) 

Xl  +  N-M) = J  1ND( NMAX  +  N-M) 

Xl+N-M) =K IND (NMAX+N-N ) 


J 
'K 


7)  NMAX,  ITEXT(IND) 

Hl//^X,12hTHE  HIGHEST  ^IS^l-^rH  HARMONICS  OF 

i  N  M  A  X 

) 

) 

1 

IPUl)  GO  TO  9 

1-1 

10)  ITEXT{^), JliKl 


A7) 


,10)  ITEXT (5),J1,K1 

/2X>A3,2X,2HJ':>I5/2X>2HK' 

E 

>12)  Y(N) 


>I5//) 
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12 

8 


21 


22 

20 


31 
C 


32 


321 


FDR 

COM 

DO 

DO 

DO 

F  VA 

IF( 

IF( 

IF( 

IF( 

IF( 

IF{ 

IF( 

IF( 

IF( 

IF( 

1F( 

IF( 

DO 

F(  I 

CON 

COM 

II- 

1L- 

MGV 

DO 

DO 

DO 

F(  I 

DO 

DO 

I3« 

I<»« 

DO 

XI- 

X2- 

X3« 

XA- 

F(I 

F(  I 

F(  I 

F{  I 

CCS 

I9« 

DO 

DO 

F{  I 

DO 

DO 

IE2 

lEJ 


MAT(10X,8HL2-N0RM«,£20.10) 

TINUE 

20  J«1#I1 

20  K«1,I2 

21  I«IFS*LF 
L{I)-F(I>J* 
ISD.EQ.l)  C 
IND.E 
IND.E 
INO.E 
IND.E 


IND 
I  NO 
IND 
IND 
IND 
IND 
IND 
22 


.E 
.E 
.E 
.E 
•  E 
.E 
.E 
I« 


>  Jf  K) 

TINUE 
TINUE 
1 

IMAX 
E  EAC 
31  1  = 
31  K» 

31  J  = 

>  J*  12 

32  J" 
32  K» 
12  +  2- 
IPU1  + 
32  I« 
F(  I,  J 
-F(I* 
F  (I>I 
-F(  I, 
*  J»  13 

>  J^K  ) 
»  I^^K 
,  IA»  I 
TINUE 
12*1 
321  I 

321  K 
,  IPUl 

322  I 
322  J 
-IPVl 
•IPUl 


If 
l» 
IF 
F( 

0. 

0. 
0. 

0. 
0. 
0. 
0. 
0. 
0. 
Q. 
Q. 
0. 
1/ 
-P 


2) 
3) 
1  ) 
1) 
1) 
2) 
2) 
2) 
3) 
3) 
3) 


S 

K) 
ALL 
ALL 
ALL 

MAX» 
MAX" 
f^AX" 
MAX« 
MAX  = 
MAX" 
MAX" 
MAX=» 
MAX" 


SINrRO{FVAL*PC*IMAXRO*NI) 

SlNT(FVAL*PC*II1AXPSIfrU/IND) 

SINT{FVAL,PC/IMAXPSI*NI^IND) 

IMAXRG 

JfXkG 

KMXRG 

IMAXPSI 

J  M  X  P  5  I 

KMXPSI 

IMAXPSI 

JMXPSI 

K  M  X  P  S  I 


IMAX 
C(I  ) 


H  K-CDLUMS  0VbR>F0R  K"IPV1  TO  IPVl+IPV 

II>  IL 

1/IPVl 

l/Il 

+2-K) =F( I> J> 12+1-K) 

1  *  I  P  U 1 

l^IPVl 

K 

J 

#K)+F( l,^, 13) 
I^*K)+F(I*  I^>  13) 
^^K)+F(I, I^, 13) 
J* 13) +F( I, J>K) 
)"X1 
"X2 
)«X3 
3)"X'» 


«II»  IL 

«1»I9 

+1#K) "0, 

»II>  IL 

"1/IPUl 

♦1 

+  J 
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322 
C 


33 


3^ 


35 


36 


F(I 

F(  I 
OUT 
IF( 
I3» 

DO 

DO 

DO 

IF( 

CON 

DO 

DO 

DO 

KK« 

IF( 

CON 

I5» 

I6» 

DO 

DO 

DO 

KK» 

IF( 

CON 

DO 

DO 

DO 

IF( 

CON 

RET 

END 


* J>IE2)-0.5*F( I, J*IE2) 

>IEJ*IPV1)-0.5*F(I^IEJ*IPV1) 

PUT  TO  TAPE82  OF  ARRAYS  CC>CS*SC*SS 

IS0.LE.3)  WRITE(82)  IMAX, J  MAX, KM  AX 

IPVl+1 

IFVl+KMAX 

33  I»II>IL 

33  J»1>JMAX 

33  K-I3,U 

IND.LE.3)  WRITE{82) 

TINUE 

3^  1= II»  IL 

3^  J«1,JMAX 

3^  K»1,KMAX 

IPVl+l-K 

WRITE (82) 


F(  I>J>K) 


INC.LE .3) 
TINUE 
IPUl+1 
IPUl  +  J.rAX 
35  1=11,  IL 
35  J  =  I5,  15 

35  K«1,KMAX 
IPVl+l-K 

IND.LE.3)  WR1TE(82) 
TINUE 

36  1  =  11,  IL 
36  J-I5,  16 
36  K»I3,  I^ 
IND.LE.3)  WRITE(62) 

URN 


F(  I,  J,KK  ) 


F(  I» J,KK  ) 


F( I, J,K) 


SUBROUTINE  SINT(  F,e,  iMAX,NI,  IND) 
C   THIS  SUBROUTINE  APPROXIMATES  IN  S  THE  DATA 
C   FOR  DBGRADU  AND  D3GRA0V  dY  POLYNOMIAL  LEAST 
C   SQUARES  METHOD. 

C0MMCN/VAR3/PI2,EP,t.'LZ,ALF 

DIMENSION    A(10,lC'),F(13),B(ll.  ) 

h.  1  1 «  M  - 1 

HS«1.0/FLl:aT(NI1) 
C   COMPUTE  THE  INNER  PRODUCTS  FOR  THE  ENTRIES 

DO  10  J  =  l,  IMAX 

DO  10  K-J,  IMAX 

NEXP-J+K-2 


OF  MATRIX  A 
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11 


10 

c 


21 

20 
C 


31 


30 


100 


SUM« 
IF(  I 

DO  1 
SL-S 
SALF 

SUM- 
CON  T 
A(  J« 
A(K> 
CDNT 
CCMP 
DO  2 
NEXP 
SUM« 
IFd 
DO  2 
IF{  I 
SL«S 
SALF 
SUM» 
CONT 
B(  J  ) 
CONT 
SOLV 
KE»1 
IA-1 
IB«1 
CALL 
Il-I 
SUM" 
IFd 
DO  3 
IFd 
IF(  I 
SL  « 
SALF 
VAL  = 
DO  3 
VAL- 
T  « 
SUM. 
CONT 
SUM« 

hKlT 

FORM 
RETU 
END 


E.3)  SL«0.0 

M»1>NI1 

*ALF 

(SALF**  f<  EXP) 

UM 
UM 

THE  RIGHT  HAND  SIDE  OF  AX-6 
Id  MAX 


E.3)  SL>=0.0 

M»1,NI1 

E.3)  EF«F(LAM+1) 

*ALF 
EF*(SALF**NEXP) 


0.0 

ND.L 
1  LA 
L+HS 
-SL* 

SUM* 
INUE 
K)=.S 
J)-S 
INUE 
UTE 

0  J" 
-J-1 
0.0 
ND.L 

1  LA 
ND.L 
L+HS 
=  SL* 
SUM  + 
INUE 
=  SUM 
INUE 

E  AX«B 

0 
0 

LEQ(A*BdMAX,NEdA>I3>DET) 
MAX-1 
0.0 
N0.LE.3)     SL=0.0 

0  I=2>NI 

ND.LE.3)     EF=F ( I ) 
ND.LE.3)     IY=1 

SL  +  HS 
.SL**ALF 
B(  IMAX  ) 

1  J  =  li  U 

VAL*SALF+B(IMAX-J  ) 
tF-VAL 
SUM+(T*T) 

INUE 

S0RT(SUM/FL0AT(NI1) ) 
E(odOO)  SUN 
AT(10X,8HST.DEV 
RN 


-/E20.10) 
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C 
C 


12 


11 
10 


15 


I'. 


20 


SU 

TH 
FO 
CO 
DI 
If 
hi 
HI 
HS 
DO 
DO 
NE 
SU 
SL 
DO 
SL 
SA 
SU 
CD 
A( 
A( 
CO 
CO 
A( 
DO 
SU 
SL 
NE 
DO 
SL 
SA 
SU 
CO 
B( 
CO 
B( 
NE 
lA 
IB 
CA 
PC 
DC 
PC 

cc 

PC 
PC 
SU 


BRDUTI 
IS  SUB 
R  PO  B 
MMON/V 
MEKSIO 
AX«  IMA 
1-Nl-l 
2-NI-2 
•1.0/F 

10  J- 

11  K» 
XP» J+K 
M»0.0 
•0.0 

12  LA 
«SL+HS 
LF'SL* 
M=SUM+ 
NTINUE 
J*K)»S 
K,  J  )-S 
STINUE 
NTINUE 
1>1)«A 

1^  J' 
M=C.0 
"0.0 
XP«J-1 

15  LA 
•SL+HS 
LF«SL* 
M  »  S  U  .i  + 
NTINUE 
J ) «SUM 
NTINUE 
1)*B(1 
»1 
•10 
-10 

LL  LEO 
( l)»-3 

20  I» 
(  I)»3  ( 
NTINUE 
(  IMAXl 
(2)«PC 
M«C.O 


NE  SINTPOCF, PC/  IMAXl, NI  ) 

ROUTINE  APPPQXIMATES  IN  S  THE  DATA 

Y  POLYNOMIAL  LEAST  SQUARES  METHOD. 

AR3/PI2, EPfOLZf ALF 

N  A{10/lC)>F(i3),B(10),PC(10) 

Xl-1 


LDAK  Nil) 
1*  IMAX 
J,  IMAX 
-2 


M=1,NI2 

*ALF 
(SALF-1.0)*(SALF-1.0)*(SALF**NEXP) 


UM 
UM 


(1/1)  +1.0 
1/  IMAX 


M«1/NI2 

*ALF 

(F(LAM  +  1)-F{M)*SALF)*(SALF-1.C)*(SALF**NEXP) 


)-F(l) 


( A/ Q,  IMAX/NE»  lA/ I3/DET) 

(  1  ) 

2  *  I  ^.  A  X 

I-l )-B(I ) 

)  •=  B  (  I  M  A  X  ) 
(2 )  +  F (NI  ) 
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31 


30 


SL  — HS 

DO  30  I«1>NI1 

SL  -SL+HS 

SALF«0.0 

IF(SL.GT. 0.001)  SALF«SL*ALF 

VAL-PC( IMAXl) 

DO  31  J-1*IMAX 

VAL»VAL»SALF+PC( IMAXl-J) 

T       -fCn-VAL 

SUM»SUM-KT*T  ) 

CONTINUE 

SUM'SORTCSUM/FLOAKNID) 

RETURN 

END 


SUB 

c 

THI 

c 

THO 

COM 

COM 

COM 

COM 

COM 

IF{ 
IF( 
IF( 

IF( 

IF( 

II- 

12- 

IMA 

DO 

DO 

XJ- 

XK- 

JP- 

KP- 

00 

IPl 

IP2 

Fl( 

Fl( 

Fl( 


ROUTINE     DaGPADS(IND) 

S    SUBROUTINE    COMPUTES    THE    COEFFICIENTS    OF    DBGRADS    FROM 

SE    OF     DBGf-  AOU    AND    03GRADV. 

MON/VAPl/M^NJ^NK^JVAC^NNDEG 

MCN/VAK3/P  l?,\:P,:iiZ,ALF 

MGN/VARB/IPU/IPV,IPU1/IPV1 

MCN/VAR21/Zi'^l(52)^Zi2(52)^Z.'^3(t>2>:>2)>ZM't('j2>52)iZM5(52,p2)/ 

Fl(13>26^26)^  F2(13/26*26) 
M0N/VAR22/IMAXR0»1MAXPSIW^MAX> 

KAXIS*JMXRO>KMXPO>JMxPSI>Ky,  XPSI 
IND.EQ.l)  S0-( 1.0/FL0AT(NI-1 ) )**ALF 

SO-C.O 


IMAX-IMAXPSI 
J  M  A  X  »  J  M  X  P  S  I 

K,''AX  =  KMXPSI 


IND.EQ.l) 
IND.EQ.l) 
IND.EQ.l) 
IND.EQ.l) 
IPUl+IPUl 
IPV+IPV+1 
Xl-IMAX+1 
1  J-l^IPUl 

1  K-1,IPV1 
J-1 

IPVl-K 
IPUl+J 
I2+2-K 

2  I«l»  IMAX 
»  I  M  A  X  1  -  I 
-IMAXl+1-I 

IP2* J  ,KP) -XJ*F1 (IP1» J 
IP2,J  ,K  )  »XJ*F1(  IPl, J 


,KP  )/FLOAT(  IPl ) 
>K  )/FLUAT{IPl) 


IP2*JP,KP)»XJ*F1(IP1,JP>KP)/FL0AT(IP1) 
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Fl( 

F2( 

f2( 

F2{ 

F2{ 

CON 

Fl{ 

Fl( 

Fl( 

Fl( 

F2( 

F2{ 

F2( 

F2( 

Tl- 

T2- 

T3- 

T^" 

T5- 

T6« 

T7- 

18- 

00 

Tl- 

T2» 

T3- 

T4- 

T5- 

T6- 

T7- 

T6« 

CON 

Fl( 

FK 

Fl( 

FK 

F2( 

F2( 

F2( 

F2( 

DO 

Tl- 

T2- 

T3» 

T<i- 

FK 

FK 

FK 

FK 

CON 

CON 

CUT 


IPZ,J  » 
IP2,J  , 
IP2,  JP> 
IP2,JP# 
TINUE 
1*J  *KP 
1*  JP,KP 
1>J  *K 
IfJPpK 
1,J  *KP 
1,  JP>KP 
1*J  >K 
1>  JP»K 
FK  IMAX 
F  K  I  .^  A  X 
FK  IMAX 
FKIMAX 
F  2  {  I  M  A  X 
F2( IMAX 
F  2  (  I  M  A  X 
F2(IMAX 
3  I-KI 
T1*S0+F 
T2*S0+F 
T3*S0+F 
TA*SO+F 
T5*S0+F 
T6*S0+F 
T7«S0+F 
T8*S0+F 
TINUE 

J 

JP 

JP 

J 

J 

JP 

JP 

I« 
-FK 
-FK 
FK 
FK 
I»J 
I»J 
I,  JP 
I>  JP 
TINUE 
TINUE 
PUT     TO 


K  )»XJ*FK  IP1> JP,K  )/FLOAT(IPl) 
KP) -XK  +  FZC IPli J  »KP  )/FLOAT(  IPl ) 
K  )»XK*F2(  IPl/J  ,K  )/FLOAT(IPl) 
KP)-XK*F2(IP1*JP»KP)/FL0AT(IP1) 
K     )  »XK»F2(IP1>JP*K     )/FLOAT(IPl) 


1* 
1^ 
1* 
K 
I, 
I, 
I, 


K 

KP 

K 

X.P 

K 

KP 

K 

KP 

*I 

>  J 
»  J 

>  J 
*  J 
KP 
K 
K 
KP 


)-0.0 

)«0.0 

)-0.0 

)-0.0 

)»0.0 

)«0.0 

)»0.0 

)«=0.0 

1>J     »K     ) 

1>J     »KP) 

l^JP^K     ) 

1*JP,KP) 

1/J     M     ) 

1»J     >KP) 

1*JP^K     ) 

If JP/KP) 

MAX 

1(IMAX1-I»J     ,K     ) 

1(IMAX1-I>J     ,KP) 

1 (  IMAXl-I*  JP»  K     ) 

1(IMAX1-I> JP>KP) 

2(IKAX1-I,J     ,K     ) 

2( IMAXl-If  J    /KP) 

2  (IMAXl-I/ JP/K     ) 

2(IMAX1-I> JP,KP ) 

)»-Tl 

)»-T2 

)«-T3 

)«-TA 

)"-T5 

)«-T6 

J^-T? 

)--T8 

MAUI 

P»K  )-F2( I>  J  ,K  ) 

PiKP)  +  F2 (  1/ J  /KP) 
/<P)-r2(I/JP/NP) 
»K     )+F2(  1/ JP/K  ) 

)-Tl*PI2 

)»T2»PI2 

)»T3*PI2 

)«=T^»PI2 


TAPE82  OF  ARKAYS  CC/CS/SC/SS 
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I3-IPV1+1 

I^«IPVl+KhAX 

IF(IND.Ea.l)  WRITE(82) 

DO  33  I-1*IMAX1 

DO  33  J»1*JMAX 

DO  33  K=I3^  I^ 

IF(IND.EO.l)  WRITE(82) 
33  CONTINUE 

DO  3^  I«=liIMAXl 

DO  3A  J»1*JMAX 

DO  3A  K=1>KMAX 

KK-IPVl+1-K 

IF(INO.EQ.l)  WRlTE(e2) 
3^  CONTINUE 

I5»IPUl+i 

I6-IPU1+JMAX 

00  35  I»1^IMAX1 

DO  35  J=I5# 16 

DO  35  K»1*KMAX 

KK«IPV1+1-K 

IF(IND.EO.l)  WRITE(82) 

35  CONTINUE 

DO  36  I=1>IMAX1 
DO  36  J«I5^  16 
DO  36  K«I3>  I^ 
IFdSD.EQ.l)  WRITE(82) 

36  CONTINUE 
RETURN 
END 


IMAX1> JMAX*KMAX 


F1(I>J>K) 


Fl( I/J^KK) 


Fl(  I* J>KK) 


F1(I>J/K) 


SUBROUTINE  SUBHJK 
C   THIS  SUBROUTINE  CO'IPUTES  DKGRADU  AND  DKGRADV  AT  (I»J>K) 
C   BY  CALLS  TO  THE  SUBROUTINE  EVAL.  AT  EACH  FIXED  SfTHESE 
C   DATA  ARE  THEN  FOURIER  ANALYZED  bY  CALLS  TO  THE  SUBROUTINE 
C   CFORf^..  THEN  FITTING  IN  S  IS  ACCOMPLISHED  BY  CALLS  TO 
C   THE  SUBROUTINE  CSINT. 

CCiMON/VARl/NIW^J^NK^NVAC^IOEG 

C0M^DN/VAP2/HS^hU>HV  i 

CO^MON/VAi:.3/PI2>  EP,  CLZf  ALF 

C0MMDN/VAR20/JMAX,Kr,AX>SO/NRES 

COMMON/ VA  R21/ Z,- 1  (52  )^Z^',2(  52)^  Fl(  52/52  )>F2{52i  52  )>F3  (52^52)* 
+  ZM3(13, 26, 26)/ZM^( 13*26,26) 

IOER-2 

IU»0 

IV«0 
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S-SO-HS 

DO  1  I-1*NI 

S-S+HS 

SL»S**ALF 

U--HU 

D  0  2  J  «  1  ^  N  J 

U»U+HU 

V--HV 

DO  2  K-1,NK 

VV  +  HV 

CALL  EVAL(SL>U^V/  I  DE  R*  X,*10D6/ DKS  /  DKU  ,  DKV  ) 

Fl( J,K)=DKU 

F2( JiK)»OKV 

IF  (I  .EQ.NRES)  F3(J^K)=DKS 

CDNTISUE 

CALL  CF0RM{ lU/IV^Fl) 

CALL  CFOPr.  (IU*IV»F2) 

M«8 

CALL  WRIT69(M,F1) 

CALL  WRITb9(M,F2) 

CONTINUE 

CALL  CFORM( IU>IV>F3) 

M«8 

CALL  CSINT(M) 

RETURN 

END 


SUBROUTINE  WR  ITe9( M, F ) 
C   THIS  SUBROUTINE  WRITES  ONTO  TAPEM  (M  =  8  OR  ""  )     FOR 
C   TE'iPDRARY  STORAGE.  IT  IS  CALIED  BY  SUBROUTINE  SUBHJK. 

CDMM0N/VAk8/IPU^IPVWPU1^IPV1 

DIMENSION  F  (52/52  )»A(p2 ) 

I1«IPU1+IPU1 

I2«IPV1+IPV 

DO  1  J-1* II 

DO  2  K»l,  12 
2A(K)-F(J/K) 
1  W  R  I  T  E  (  M  )  A 

RETURN 

END 
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SUBROUTINE  PEAD89(M, J,F1,F2) 
C   THIS  SUBROUTINE  READS  TAPEM  (M^S  OR  9)  AS  THE  DATA 
C   STORED  ON  THEM  ARE  NEEDED  BY  THE  SUBROUTINE  CSINT. 

COMMON/ VA  R  1 /M  #  N  J,  NK,-N  VAC*  IDEG 

C0MMCN/VAR3/IPUWPV>  IPJl*  IPVl 

DIMENSION  Fl( 52^52)* F2(52> 52 ),A{ 52) 

Il-IPUl+IPUl 

I2-IPV1+IPV 

I3-IPU1+IPU 

NI-NI 

REWIND  M  :'  :  ,  ,.   .-   *  .  ^     V  ,  •  .  .  :  - 

IF(J.EO.l)  GO  TO  1 

Jl-J-1 

DO  2  JJ-1^  Jl  ....:.  /;  -,  ■ 

2  READ(M) 

1  CONTINUE 

DO  3  I«1WU 

IF(I.GT.l)  GO  TO  31 

GO  TO  32 

DO  33  JJ«1*  13 

READ(M) 

CONTINUE 

READ(M)  A 

DO  4  K»1,I2 
A  F1(I,K)=A(K) 

DO  5  JJ=1*  13 

5  READ(M) 
READ{M)  A 
DO  6  K»1,I2 

6  F2( I,K)»A(K) 

3  CONTINUE 

RETURN  -  .  .      , 

END 


31 
33 
32 


SUBROUTINE  CSINT(M) 
C   THIS  SUBROUTINE  HAS  THREE  FUNCTIONS: 

C     1.  IT  PCLYNO.^IAL  LEAST-SQUARES  FITS  DATA  FROM  THE 
C         SUbRCUTiSE  SU3HJK. 

C     2.  IT  COMPUTES  THE  COEFFICIENTS  OF  UKGRADS  FRON 
C         THOSE  OP  OKGRADU  AND  DKGRACV, 
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C     3.  IT  OUTPUTS  THE  FINAL  COEFFICIENTS  OF  DK GR ADU> DKGR ADV * 
C        AND  DKGRADS  TO  TAPE83  FOR  USE  BY  THE  LINE/ORblT  CODE 
C        LIN0R3.  TO  DO  THIS  IT  CALLS  SUBPOLTINE  WRITa3. 

COMMON/VARl/NI,NJ,NKiNVAC»IDEG 

C0MMCN/VAR2/HS»HU/HV 

C0M10N/VAk3/P  12, EP» OLZ» ALF 

C0MM0S/VAkS/IPU,IPV,IPU1,IPV1 

C0MMCN/VAP20/JfiAX»KMAX,S0>NRES 

COrt.iON/VAR21/Al(^2)>A2(52)>Fl(52,52),F2(52,52),F3(52>52), 
+  C0Ul(13,2o>26)*C0Vl (13>26>26) 

DIMENSION    COUdOi  20  ),CJV(iO>20)K(10^1G),Cl{10>10), 
+  F<.(26/20) 


101 


100 


c 

cou  « 

GOV    « 

1-IPU 

2»IPV 

1«NI 

RES«S 

LO-SR 

OMPUT 

0    ICO 

0  100 

EXPi:J 

UM1«0 

1  -S 

0  101 

1  «s 

L  »S 
UM1  =  S 
(  J»K) 
(  <  ,  J  ) 
CNTIN 
IND=1 

1  «1 

2  «J 

1  -I 

2  «K 
QNTIN 
0  2  J 
F  (  N  I  N 
F  (  N  I  N 
ALL  R 
G  3  K 
0  ^  I 
1(  I  )  - 
2(1)- 
ALL  L 
ALL  L 
CNTIN 
0  30 
0  30  ■ 


10 
10 
10 
1  + 
1  + 

0  + 
ES 
E 

J 

K 

+  K 

.0 

0- 

L 

1  + 
1* 

u.i 
=  s 

-S 

UE 


IPUl 
IPV 

(NRES-1.)*HS 

♦*ALF 

THE  RIGHT  HAND  SIDE  OF  CX«tJ 

-1,  IDEG 

=  J,  IDEG 

-2 

HS 

AM=1, Nl 

HS 

♦  ALF 

1+SL**NEXP 

U.'il 

urii 


MAX 

PVl 

MAX 

JE 

'Jl 

D.E 

D.  £ 

EAj 

'I, 

«1, 

Fl( 

F2( 

HS( 

HS( 

UE 

IJ- 

IK- 


-KMAX 
-1+KMAX-l+l 


>  J 

C  . 
G  . 
c9 
^2 
Nl 
1/ 
I, 


2 

1)  JU«J 

2)  JU='J-J1  +  JMAX  +  1 
( ."^  *  J  /  F  1  >  F  2  ) 


K+Kl  ) 

K  +  Kl) 

IDEGiM,  A1,C0U) 

IDEG>N1, A2,C0V) 


If  IDEG 
1/IDEG 
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30  CKIJ,  IK)-C(IJ>IK) 

CALL  LEQ(C1>CCU*I0EG»K2*IC*IC0U*DET) 
DO  31  IJ»1,IDEG 
00  31  IK»lfIOEG 

31  Cl( IJ*IK)«C (IJ/IK) 

CALL  LEQ(C1,CGV*IDEG*K2,IC>ICDV»CET) 

DO  5  K-1^K2 

DO  5  I-l^ lOEG 

COUK I* JJ*K)-Cau(  I»K) 

5  CDVK  I,  JU>K)»CCV(  It  K) 
2  CONTINUE 

IF(NIND.EQ.2)  GO  TO  6 

MND«2 

Jl   «IPU1+1 

J2   »IPU1+JMAX 

GO  TO  1 

6  CONTINUE 
JFl-1 
JF2«JMAX 
KF1«IPV1-KMAX 

KF2-KMAX-1+KMAX-1+1  .  •      :    r 

DO    7     J-JF1,JF2 

DO    7     K=1>KF2 

F^( J*K)«F3( J/K+KFl) 

F'.CJ  +  JFZ^Kj^FBCJ  +  IPUl^K+KFlJ 

7  CONTINUE 
J2-JMAX+JMAX 
K2«KMAX+KMAX 
DO    10    I-1,IDE6 
DO    10    J»l>  J2 
DO    10    K»1,KMAX 

C0U1(I/J,K2  +  1-K)=CUU1(  I,  J/K2-K) 
C0V1(I,J,K2+1-K)=C0V1(I^J,K2-K) 
IF(I.NE.l)     GO     TO    10 
F^(J,K2  +  1-K)-F',(  JiK2-K) 
10    CONTINUE 

00    11     I-l/IOEG 

DO    11     J-1, JMAX 

DO     11     K  =  1,N,MAX 

Xl«CGUl(I,J,K2+l-K)+CJUl(IfJ,K) 

X2-CCUl(I,JMiX+J,K2+l-K)-CLUl(I,JMAX+J,K) 

X3»CCU1(1>J'<AX+J,K2+1-K)  +  CCU1(1^JMAX  +  J,K) 

X't-CGU1(1>J,K)-CJU1(I>J>K2  +  1-K) 

COUK  1/ J,K2  +  1-K)  «X1 

C0U1{I,J,K)  .X2 

COUl (  I> JMAX+J,K )  =X3 

C0U1(I,JMAX  +  J,K2>-1-K)»X^ 

X1-CCV1(I,J,K2+1-K)+C0V1(I>J>K) 

X2"C0V1(I,JMAX+J,K2+1-K)-CGV1(I/JMAX+J,K) 

X3»C0V1( I*JMAX  +  J,K2  +  1-K)+CQV1 (I j JMAX* J,  K) 

X<.-C0V1(I>J,K)-C0V1(I^J,K2  +  1-K) 

COVK  I»  J*K2  +  1-K)  -XI 
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=  0.0 
'0.0 


AX  ) 


»0.  t 

-0.5 
=  0.5 
»0.  5 


C0V1(I#J>K)  "XZ 

COVK  I*  J.MAX+J,K)  «X3 

C0Vl(I>JMAX  +  J,K2  +  l-K)«X'i 
IF(I.NE.i)     GO     TO     11 
X1»F4(J*K2  +  1-K)+F4(  JM) 
Xa-FACJ'iAX  +  JfKa+l-KJ-F^CJ 
Xa-F'fCJMAX+JtKZ  +  l-Kj+F^CJ 
X<»»F^{  J*K)-F^  (  J^K2+1-K) 
F<»(Jf  K2  +  1-K)  »X1 

F^(J,K)  «X2 

F^(  jr'AX  +  j^K )  -xa 

F^( JNAX+J>K2+1-K) «X^ 

11  CONTINUE 

DO    12     I-l> IDEG 
DO    12    K-1, K2 
FA(JMAX+1*K)«0.0 
COUK  I,  JMAX  +  l^K) 

12  CDVK  I>  J-.AX  +  liK) 
DO  13  I»l,  IDEG 
DO  13  J=lfJMAX 
COUK  I*  J>KMAX  +  1  ) 
COUK  I*  JMAX+J,KM 
COVl (  I^  J>KKAX  +  1) 
COVKI>JMAX  +  J,KM 
IF(I.NE.l)  6G  TO 
F<i(  J,KMAX  +  1  )  =0,5 
F^( JhAX+J>KMAX)= 

13  CONTINUE 
C       OUTPUT    TO    TAPE83 

CALL    WRIT83(C0U1 
CALL    WRIT&3(C0V1 
I0EG1»IDEG+1 
DO    20    J«1,JMAX 
DO    20    K«1>K,MAX 
XJ»J-1 
XK-KMAX-K 
JP-JMAX+J 
KP-K2+1-K 
DO    21     I«l» IDEG 
IPl-IDE&l-I 
IP2«  IDEGl+1-I 
COUK IP2^ J>nP  )     «XJ«CDU1(I 
COUl  (  IP2,  J,K)        =X  J^CGUK  I 
CGUK  iP2,  JP>  ■<  P)  =  XJ*CCUi(  I 
COUK IP2^ JP>K )     »XJ*CDUl(i 
COVK  IP2^  J>KP  )     -XK*CGVK1 
C0V1(IP2,J^K)        =XK*C0V1(I 
C0VKIP2,J?>KP)-XK*C0VKI 
COVK IP2f JP>K )     -XK*CuVl(I 
21    CONTINUE 

COUKl^J^KP)     «0.0 
COUK!*  JP^KP)  "0.0 


.i  A  X  +  J  i  K  ) 
MAX+  J, K) 


AX  ) 
13 
*^  F  -i  (  J  ^  K  M  A 
0.5*F^( JM 


»  1  D  E  G  ^  J  M  A 
,  IDEG, JMA 


«CDU1(  I>J>KMAX  +  1) 
«C0U1(  I* JMAX  +  J,KMAX ) 
*C0V1(  I^ J^K^AX  +  l) 
♦CUVK  K  JiMAX  +  J^KMAX  ) 

X+1) 

AX+J^KMAX) 


X  ^  K  M  A  X  ) 
X  ,  K  f^  A  X  ) 


Plf J,K  P) /FLUAT ( IPl  ) 
Pl> J  ^K ) /FLOAT (  IPl ) 
Pl>  JP^  KP) /FLOAT (IPl) 
Plf JP^K)/FLUAT ( IPl ) 
P1>J>KP)/FL0AT ( IPl ) 
PI*  J>K  )  /FLOAK  IPl) 
P1»JP/KP)  /FLOAT{ IPl ) 
P1,JP*K)/FL0AT (IPl) 
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CDU1(1*J»K)        -0.0 

C0U1(1*JP>K)     »0.0 

C0V1(1,J*KP)     »0.0 

COV1(1,JP*KP)»0.0 

COVKl^J^K)        «0.0 

C0V1{1,JP>K)     «0.0 

T1«C0U1 ( IDEG1»J^K) 

T2-C0U1(  IDEG1>J*KP) 

T3«CCU1(  IDEGlf JP^K) 

T^«CCU1 (IDEGl/ JP,KP) 

T5-C0V1( IDEGl^J^K ) 

T6»CDV1( IDEGl^ J,KP) 

T7=CDV1( IDEG1> JP,K) 

T8-C0V1( IDEGl^JPiKP ) 

DO  ZZ     I«l*  IDEG 

Tl-Tl*SLO+CQUl(IDEGi-I/ J,K) 

T2»T2*SL0+C0U1(ICEG1-I,J,KP) 

T3='T3*SL0  +  C0U1(IDEG1-I,JP,K) 

T^«T<r*SLO  +  COUl(IDtGl-I*JP,KP) 

T5"T5*5L0+C0V1(IDEG1-I,J,K) 

T6«T6*SL0+CCV1(IDEG1-I*J^KP) 

T7«T7*SL0+CaVl(IDGGl-I,JP>K) 

T8»Te*SL0  +  C0Vl(IDEGl-I,JP,i<P) 
ZZ    CONTINUE 

CGU1(1,J>K)   =-Tl 

C0U1(1>J>KP)     »-T2 

CGU1(1,JP,K)     --T3 

CGUK  1,  JP,KP)  «-T^ 

CGV1(1>J,K)        --T5 

COVKl^JiKP)     »-T6 

CCVKl^JP^K)     «-T7 

COVl ( li  JP/KP ) «-T8 

DO    23    1=1, IDEGl 

T1«-CGU1(  1> JP,K)-CDV1{ I* J,K) 

T2«-C0Ul(I>JP,KP)+Cavi{I>J>KP) 

T3»    CGU1(I,J,KP)-Cavi(  I>  JP.,KP  ) 

T^»    COUl  {  I,  J,  K)  +  C0V1  (  IOP,K) 

COUl ( I>  J>KP )       »T1*PI2 

CGU1(I,J>K)  •'T2  +  PI2 

COUl  (  I,  JP,K)        =«T3*PI2 

CCUK I, JP,KP)     =T^*PI2 
23    CONTINUE 

CGUl  (1,  J     ,KP)  -CDUK  1,  J 

CCU1(1,J     ,K    )=CGU1{1/J 

C0U1(1>J?,K     )  »C0U1(1>  JP,K     )  +  F-^(JP,K     ) 

C0U1(1,J?,KP)»CGU1(1,JP,KP)+F4(JP,KP) 
20  CONTINUE 
C   OUTPLT  TO  TAPE83 

CALL  WRIT33(C0U1,IDEG1,JKAX,KKAX) 

RETURN 

END 


,KP)+F^( J 
,K  )  +F^( J 


>KP) 
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SUBROUTINE  WRIT83(F,IMAX>JMAX,KMAX) 

THIS  SUBROUTINE  -^  R  I  T  E  5  ONTO  TAPE83. 

DIMENSION  F(13*26,26) 

I3«Kr<AX+l 

I<.»Kf,AX  +  KMAX 

WRITE(83)  IMAX> JMAX,KMAX 

DO  1  I"l^  IMAX 

DO  1  J-1>JMAX 

DO  1  K«I3> I^ 

WRITE(e3)  F(IiJ>K) 

00  2  I-I/IMAX 

DC  2  J=l>  JMAX 

DO  2  K=1^K,MAX 

KK«Kh^AX  +  l-r<. 

WRITE(83)  F(I^J^KK) 

I5=J^AX+1 

I6«JrAX+J^AX 

DO  3  I«1>IMAX 

DO  3  J=I5^ 16 

DO  3  K«1/KMAX 

KK»KMAX4-1-K 

WRITE(d3)  FdiJ/KK) 

DO  ^     1=1/ I  MAX 

DO  ^  J=I5>  16 

DO  ^  K»I3>  I^ 

WRITE(B3)  FdiJ/K) 

RETURN 

END 


SUBROUTINE  LH5(N/IDEG/N1,A,3) 
C   THIS  S-J-3RGUT1NE  COMPUTES  THE  INNER  PPQDUCTS  OF  Y 
C   WITH  POLYNOMIALS  IN  S  FOr  THE  LINEAR  SYSTEM:AX=Y. 
C   IT  IS  CALLED  BY  SUBROUTINE  CSINT. 

CGMK0N/VAfc2/HS/HU,HV 

C0MH0N/\/AR3/PI2*EP/QLZ>ALF 

C0MfiDN/VAf20/JMAX,KMAX,S0»NPES 

DIMENSION  A{52) /B (1C>20) 

DC  1  J»li IDEG 

NEXP» J-1 

SUM  =0.0 
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Sl 
DO 
SI 
SL 


»SO-HS 

!  LA,1«1,N1 
"Sl+HS 
«S1»*ALF 


2  SU.i«SU1  +  A(LAM  )*(SL**NEXP) 
1  B( J»N)»SUM 

RETURN 

END 


C 
C 
C 


5 
11 

10 
21 


SU3R 
FAST 
INPJ 
REPL 
COMP 
DIME 
IF  ( 
NS  « 
NR  = 
NQ  - 
SET 
PI»3 
DT  « 
IF(  ( 
ANG 
DO  5 
CN(  J 
SN(  J 
ANG 
DO  1 
IF  ( 
CGNT 
f.t  D  »  N 
NS  » 
Nk  « 
IQ  « 
ID  = 
DO  2 
DO  2 
L  « 
LP^L 
M»IO 
W  -F 
IF(N 
L-LP 
DO  2 


OUTI 

FQU 
T  AR 
ACED 
LEX 
NSIO 
N.LT 

1 

2 

H 
THE 
.I'tl 

(PI 
SN(1 
"  0. 


NE  FFOPM(N* FiX>CN,SN) 

RIEft  TRANSFORM 

RAY  F  WITH  REAL  AND  IMAGINARY 

BY  ITS  FOURIER  TRANSFORM 
F  ( 1  )  >  X  (  1  )  ,  W 

N  C  ^^  {  1 )  *  S  N  ( 1 ) 
.2)  RETURN 


SINES  AND  COSINES 

592b^358«i79 

+PI)/FLGAT(N) 

)  .EQ. 0.  )  . AND.  (SN(2 ) 


PARTS     IN    ALTERNATE    CELLS 


EQ. SIN(DT)  )  )    GO     TO    11 


)     « 
)«-S 

=     AN 

0  K 

^'0^{ 

1  SUE 
G/K 

NS* 
K 
0 
0 

2  I 
A     J 
IC+J 
+  S0 


1W>1 
COS ( ANG) 
IN (ANG) 
G  +  DT 

'     NRf  N 

NO^  K  )  .ECO) 


GO    TO    21 


l^NS 
1/NO 


(L)+F(LP)*CMPLX{CN(M*1),SN(M+1)) 
R.EQ.2)    GO    TO    2A 

6    K«3iNR 


-117- 


L 

M 
I 

W 
X 

I 
I 
I 
c 

N 
I 
D 
F 
R 
61  0 
I 
C 


26 
2^ 


22 


32 


7<t 


»  L+ND 

-  K+ID 
F  (M.GE.N)  M  ■  M-N 

«  W  +  F(L)*CM?LX(CN(M.+  1)>SN(M  +  1)) 
(lO+J)  -  W 
D  •  ID+ND 

F(IC.GE.S)  ia»IQ-N 

0NTINU5 

G  -  ND 

F  (ND.GT.l)  GO  TO  61 

0  32  K  «  1>N 

(K)  «  X{K  ) 

ETUfeN 

0  60  K  «  URfH 

F  (MOD(NQ>K)  .EO.O)  GO  TO  71 
60  CONTINUE 
71  ND«NC/K 
NS  ■  KS*K 
NR  «  K 
10  -  0 
ID  «  0 

DO  72  I  =  1>NS 
DO  7 if  J  «  1>ND 
L  -  IQ+J 
L  P  «  L  +  N  D 
M«=ID 

w«x(L)  +  x(LP)♦c^'PLx(c^(rl  +  l)>SN(^ul)) 

IF(NR.E0.2)  GO  TO  7^, 
L»LP 

DO    76    K»3>MR 
L     »     L+ND 
M     «     h  +  I  D 

IF  (.-.GE.S)  M  =  r.-N 
76  W  «  W  +  X(L  )*CMPLX(CN  (N  +  1  )>  SN(r.  +  l  )  ) 

D+J)  »  ^ 

«  ID+ND 

I  C  +  N  Q 

IQ.GE.N)  lO-IO-N 

TINUE 

=  ND 

(ND.GT.l  )  GO  TO  11 

URN 


(  I 
D 

Q- 
F( 
72  CON 

w 

F 
ET 

ND 
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SUBROUTINE  PO L Y ( A , I  I  , S > VAL  ) 
C   THIS  SUBROUTINE  COMPUTES  A  POLYNOMIAL  IN  S  OF 
C   DEGREE  II  BY  THE  HORNER  /IETHOD. 

DIHENSIGN  Ad) 

VAL-A( II+l) 

IF(II.Eu.O)  RETURN 

DO  10  I-l>  II 

VAL»VAL*S+A(II+1-I) 
10  CONTINUE 

RETURN 

END 


SUBROUTINE  F0URIER(X,Y*Ka>A,CS2,SN) 
C   THIS  SUBROUTINE  COMPUTES  A  ONE  DIMENSIONAL  FOURIER 
C   SERIES  3Y  A  METHjD  OF  J.M.UATT, 

C   IN  "A  NOTE  ON  THE  EVALUATION  UF  TRIGONOMETRIC  SERIES'* 
C   COMPUTER  JOURNAL*  VOL. l>«f^>  (1959)  . 

DIMENSION  X(l  ),Y( 1) 

TRX1«0.0 

TRY1«0.0 

TRX  »X(K0+1) 

TRY  «Y(kO+1 ) 

DO  10  K«liK0 

TX   «TRX*CS2-TRX1+X(K0+1-K) 

TY   'TRY  +  CS2-TRY1  +  Y (KO+l-K)   •,  .  •    ' 

TRXl-TRX 

TRY1»=TRY 

TRX  -TX 

TRY  =TY 
10  CONTINUE 

AX   =TKX-TRX1*0.5*CS2 

A    «TRY1*SN 

A    «^A  +  AX 

RETURN 

END 


-119- 


C 
C 


10 


11 


12 


SUB 

ThI 

IN 

COM 

COM 

COM 

COM 

COM 

COM 

DIM 

DTU 

DTV 

DO 

DO 

CMA 

CAL 

DO 

CAU 

CON 

DO 

DO 

CMA 

CAL 

00 

CAU 

COS 

CAU 

DO 

CAU 

DO 

CAU 

DO 

DO 

Al- 

A2  = 

A3 

A^« 

ANG 

F  (  J 

F  (  I 

ANG 

F  (J 

F  (  I 

COS 

DO 

F(l 

F  (  I 

CON 

F(l 

RET 

END 


ROUTIN 
S  SU3R 
U  AND 
PLEX  C 
MON/VA 
M  0  f^  /  V  A 
MON/VA 
MON/VA 
MON/VA 
ENS  ION 
"0.5/F 
-O.^/F 
3  J-1^ 
^  K«l> 
T(K)-C 
L  FFOR 

5  K=l, 
X( J*^  ) 
TINUE 

6  K«=l/ 

7  J=l/ 
T(  J)-C 
L  FFOk 

8  J=l, 
X( J^K) 
TINUE 

X  ( ru  + 1 

9  K'l, 
X(NJ+1 
IC  J'l 
X  (  J  ,  N  K 
11  J=l 

11  K-1 
REAKC 
REAKC 

-  A  I  M  A  G 

-  A  I  M  A  G 
»PI2*( 

>  IPV  +  K 
PUl+J , 
-PI2*( 

>  I  P  V  ♦  2 
P  U 1  +  J  f 
TINUE 

12  K»l 
»K  )  «0. 
PU1+1# 
TINUE 

,  IPVl  ) 
URN 


E  CF 
OUTI 

V  0 
MAT/ 
R  1  /  N 
R3/P 
R6/C 
R7/C 
k8/I 

F(5 
LOAT 

loat 

NJ 
NK 

M  ?  L  X 
^  {  NK 

NK 
"CMA 


ORM{ 
NE  S 
F  DA 
CAUX 
I>NJ 
I2>E 
NU(  5 
MAT  ( 
PV,  I 
2,52 
(NJ  ) 
(NK  ) 


lU/ 
UP£ 
TA. 

>X 

,  NK 
P/C 
0), 
51  ) 
PV/ 
) 


IV/F) 

RVISES  THE  FOURIER  ANALYSIS 


>NVAC> IDEG 

LZ, ALF 

SNU(50)/CNV( 50), SNV ( 50) 

,CAUX( 51, 51 ),X(51) 

IPUl, IPVl 


(  F  (  J  ,  K  )  ,  0  .  0  ) 

,Ct<AT,X,CNV,5NV) 

T(K)  /FLDAT(NK) 


NK 

NJ 

AUX( J,K  ) 

M  (NJ  ,Cr  AT,X,Ci^U,SNU) 

NJ 

=CMAT( J)/FLDAT(NJ) 


,  NK  + 
NK 

,K)  = 
,NJ 
+  1)  = 
,  IPU 
,  IFV 
AuX( 
A  U  X  ( 
(CAU 
(CAU 
(  J-1 
)=A1 
IPV  + 
(  J-1 
-K  )  = 


1)  ^  C  A  U  X ( 1 , 1  ) 


CAUX (1, K) 


CAU 

1 

1 

J  *K- 

J,N 
X(N 
X(N 
.  )* 

♦  CG 
K)* 
.  )  * 
A2'!' 


X  (  J  ,  1  ) 


)+C AUX(NJ+2-J,NK+2-K ) ) 

K  +  2-K)  +  CAUX(r;J+2-J,K)  ) 

J+2-J,NK+2-K)-CAUX(J,K)) 

J+2-J,K)-CAUX(J,NK+2-K)) 

IU*DTU+(K-1. ) *IV*DTV) 

S(ANG)-A3»SIN(ANG) 

AS^'CDS  (  ANG)  +A1^SIN  (ANG) 

I'J*CTj-(K-i.  )  *iV*DTV  ) 

C3S  (  ANo  )-A'^vS  II.  (  ANG  ) 

)  =  '»t*C  JS  (  ANG)  +A2*SIN  (  ANG) 


,  IPV 

0 

K  )  "0  .0 

■=0.5*F(1,  IPVl) 
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SUBROUTINE  FSUB(CC»CS*SC/SS>II>JJ*KK,S*T»I1,J1*K1) 
C   THIS  SUBROUTINE  EVALUATES  A  SEPIES  CONSISTING  OF 
C   POLYNOMIALS  IN  S#SINE5  AND  COSINES  IN  U  AND  V. 

C0MiCN/VARll/CSU2>SINU»CSV2*SlNVf P( 10)^0(10)* PP{10)iOO(10)* 
+  PPP(IO) 

DIMENSION  CC(I1*J1*K1),CS(I1*J1*K1),SC(I1>J1>K1),SS(I1/J1>K1) 

IIl-II-l 

JJl-JJ-1 

KK1«KK-1 

DO  10  I'lfll 

DO  20  J«1,JJ 

DO  30  K=1/KK 

P(K)«CC (  I* J^K  ) 

Q(K)»CS (I/J/K  ) 
30  CONTINUE 

CALL  F0URIER(P>0,KK1,PP(J),CSV2,SINV) 

DO  ^0  K=1,KK 

P(K) »SC (  I*J>K  ) 

Q(K)»SS(I>J>K) 
40  CONTINUE 

CALL  P0URIER(P*C>KK1,UQ(J),CSV2>SINV) 
20  CONTINUE 

CALL  F0URIER(PP^QQ*JJ1,PPP(I),CSU2>SINL) 
10  CONTINUE 

CALL  PQLY(PPP^ II1,S> VAL) 

T-VAL 

RETURN 

END 


SUBROUTINE  DFSUB(S/XS>XU/XV*CC>CSfSC»SS»II»JJ*KK^ll>Jl>Kl) 
C   THIS  SUBROUTINE  EVALUATES  THE  DERIVATIVE  OF  A  SERIES 
C   CONSISTING  OF  PCLYNOriALS  IN  S^SINbS  AND  COSINES  IN  U  AND  V. 

COMMON/VAR3/PI2>EP*(jLZ>ALF 

C0MM0N/VAR11/CSL2>SINU>CSV2/SINV*P(10),Q(10),PP{ 10)^QQ{10) ^ 
+  PPP( 10) 

DIMENSION  CC(Il*Jl»Kl)^CS(Il/Jl>Kl)^SC(IliJl/Kl),SS(Il,Jl>Kl) 

IIl-II-l 

JJl'JJ-1 

KKl-KK-1 

DO  1  I»1^II 
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11 


12 


10 


DO  2 
DO  3 
P(K) 
Q(K) 
CALL 
PP(  J 
DO  ^ 
P(K) 
0(K) 
CALL 
QQ(J 
CONT 
CALL 
CONT 
CALL 
XU«V 
DO  5 
DO  6 
DO  7 
P(K) 
Q(K) 
CALL 
DO  8 
P(K) 
C(K) 
CALL 
CONT 
CALL 
CONT 
CALL 

xv»v 

DO  ^ 
DO  1 
DO  1 
P(K) 
Q(K) 
CALL 
DO  1 
P(K) 
C(K) 
CALL 
CONT 
CALL 
PPP{ 
CONT 
II2« 
CALL 
XS«V 
RETU 
END 


J»l> JJ 

K«l*  KK 
«SC(  I* J*K 
»SS{I> J*K 

FOURIERC 
)-FLOAT(J 

K«1>KK 
«-CC (1,J* 
--CS(  I^ J^ 

FOURIER( 
)-FLOAT(J 
INUE 

FOURIER( 
INUE 

POLY(?PP 
AL*PI2 

I«=1*II 

J=l>  JJ 

K»l*  KK 
«  FLOAT(K 
'-FLOATCK 

FOURIER( 

K>=1,KK 
=  FLDAT(K 
=-FLOAT(K 

FOURIER{ 
INUE 

FOURIER( 
INUE 

PGLY( PPP 
iL*PI2 

I«l»  III 

0  J  =  l,  JJ 

1  K=1>KK 
»CC(  I  +  1*J 
=CS(I+1>J 

FOUR  IE R( 

2  K* 1, KK 
-SC(  1*1, J 
-SS (  I+1>J 

FOUkIER( 
INUE 

FOURIER ( 
I ) -FLCAT ( 
INUE 
IIl-l 

PCLY( PPP 
AL 
RN 


) 

P»Q»KK1, PP( J)*CSV2*SINV) 
-1)*PP( J  ) 

K) 
K) 

P*0*KK1,QG(J )/CSV2>SINV) 
-1)*0Q(J ) 

PP,3C^JJ1^PPP{I)>CSU2*SINU) 

,  II1,S> VAL) 


-1  )*CS  (1^  Jj»K) 
-1)*CC(  1>  J^K) 
P/C^KK1,PP( J )*CSV2*SINV) 

-1)*SS (  I> J^K  ) 
-1)  *SCl  1/ JiK) 

P»C,KK1/Q0( J )>CSV2>SINV) 

PP^uQ,JJl>PPP{I),CSU2>SINU) 
f  IIl/S* VAL) 


P>C, KKl^  PP( J)/CSV2>SINV) 

»K) 

/K) 

p,  Q,  KKii  qcj(  J)  ,csv^^s  inv) 

PP*wC,JJl>PPP(I)^CSU2>SINU) 
1)»PPP{  1) 


,  II2^S> VAL) 
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SUBRO 

c 

THIS 

c 

THE  P 

COMMD 

COMMO 

COMMO 
t- 

COMMO 
t- 

COMMO 

COMMO 

COMMO 
COMMC 

K 

COMMO 

DIMEN 

DIMEN 

UPl 

VPl 

UP2 

VP2 

UP3 

VPS 

cu 
su 
cv 

sv 

CU2 

SU2 

CU3 

SU3 

CU3V 

SU3V 

CUV 

SUV 

CU2V2 

SU2V2 

CU3V3 

SU3V3 

CSU2 

SINJ 

CSV2 

SINV 

RAD 


UTINE  EVAL(S>U..  V*  IDER*XMOOB>DKS*DKU*DKV) 

SUBROUTINE  EVALUATES  DKGR AD S* DK GK ADU> DKGR AD V  AT 

OINT  (S,U*V) 

N/VAR3/PI2>EP>  OLZ*  ALF 

N/VARll/CSIj2>SINU/CSV2,SINV>P(10),Q(10)*PP(10)fQQ(lG)* 

PPP( 10) 
N/VAR13/IR*JR,KR,IR0*JR0>KR0>RCC(6,10^10),RCS(6,10»iG), 

RSC{6,10*10)>RSS{b>10»10) 
N/VAR1^/IB>J3/KB»I3U*JBU>K3U»CC2(6,1C,10),CS2{6,10,10), 

SC2(6,10,10)*SS2(6>10^10) 
N/VAR15/IC*JC>KC>I3V^JBV>K6V>CC3{fc>10*10),CS3(6,10,10), 

SC3{6ilO*10)/SS3(6^10,10) 
N/VAR16/IA,JA>KA*IBS>J3S,KBS,CC1{6»10,10),CS1(6,10,10), 

SC1(6^ 10>10)^SS1(6, 10,10) 
N/VAP17/KRAX,KZAX,RAC(10)^RAS(10),ZAC(10)*ZAS(10) 
N/VAkl8/DEL0>DELl,DEL2,DEL3>DEL10,DEL20,DEL30>DEL22> 

DEL33 
N/VARig/W^X.-U,  VK 

SIGN    TCC(6,iO,10),TCS(OilO,10),TSC(6,10,10)>TSS(6,10,10) 
SION    TXC (26),TXS(26) 

«PI2*U 

=PI2*V 

=2.0*UP1 

-2.0*VP1 

»3.C*bPl  ■      , 

=3.0* VPl 

«C0S(UP1) 

-SIN(U?1)     ..         ,  . 

=C0S(VP1) 

« SIN (VPl) 

=  CU*CU-Slj«SU 

=2.*CU»SU 

-CU2*CU-SU2*SU 

=SU2*CU+CU2*SU 

=«CU3*CV+SU3*SV 

»SU3*CV-CU3*SV 

»CU*CV+Sb*SV 

»SU*CV-CL'*SV 

-cuv*cuv-sav*suv 

-2.*CUV*SUV 

=  CU2V2*CUV-SLI2V2*SUV 

■SU2V2*CUV  +  C1I2V^*SUV 

-2.0*CU 

=  SU 

■2.0*CV 

■SV 

-l.-DEL0*CV+DEL10*CU+DEL20*CU2+DEL30*CU3-DEL3*CU3V 
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C 

c 
c 
c 
c 
c 
c 
c 


RADU 

RADV 
RADUU 

RADUV 

RADVV 

AN 

PN 

RAOl 

RAO 

RAD2V 

RAOV 

RADU 

RAD3 

RAD't 

RADVV 

RADUU 

RADUV 

ROB 

ZGB 

ROBU 

Z03U 

ROBV 

ZCBV 

ROBUU 

ZOBUU 

ROBUV 

ZGBUV 

ROBVV 

ZOBVV 

CALL 

CALL 

CALL 

CALL 

CALL 

DO  1 
DO  1 
DG  1 
TCC  (  I 
TCS(  I 
TSC(  1 
TSSd 
CALL 

RGSU' 
ROUU- 
RCUV- 
IROl' 
DO  2 


«P 

C  — 

c 


-3 

-1 

«R 
«  P 

S  •» 
■  — 

=  P 
»  ( 

«R 
*9 

m  — 

=  R 

=  R 

«R 

«R 

=  R 

»R 

«R 

=  R 

«R 

»R 

=  R 

«R 

FSU 

FSU 

FSU 

FSU 

DF5 

1-1 
J  =  l 
K«l 
,J> 
>  J> 
,J, 
fU 
DFS 

PQS 
RGU 
RDU 
IRQ 
i-1 


DEL22 
PI2*( 

2.*DE 
I2*(D 
PI2*P 
U3V  +  ^ 
PI2*P 
PI2*P 

0 
1.0/A 
,0  +  AN 
ADl** 
I2*Dc 
R  A  D  »  R 
3.0*P 
I2*PI 
1.0  +  A 
AD*PA 
.  0  *  R  A 
3.0*R 
AD*CU 
A  D  ♦  S  U 
ADU*C 
A['U*S 
AGV*C 
ADV*S 
ADUU* 
ADUU* 
ADUV* 
ADUV* 
ADVV* 
ADVV* 
B(CC1 
B(CC2 
3(CC3 
B(RCC 
UB(  S, 

IR 
,  IRO 
,  J5Q 
f  KRO 
K  )  .  F 
^  )  »  F 
.\  )  »-F 
K  )«-F 
UB(S> 

IR 
U*PI2 
U»PI2 
V»PI2 
-1 
f  IROl 


EL3«CU3V 
N)*RAD2V*RAD2V/RAD1+RAD3 
D4/kADl 

DVV 

ADVV 

+DtLl*CV-DEL2*CUV 

+DEL1*5 V+DEL2*SUV 

U-PI2'P.AD*SU+PI2*DEL2*SUV 

U+PI2*RAD*CU+PI2*DEL2*CUV 

U-P12*DtLl*SV-PI2*bEL2*SUV 

U+PI2*DEL1*CV-PI2*DEL2«CUV 

CU-2.*P12*RA0U«SU-P12*PI2^RAD*CU+PI2*PI2*DEL2*CUV 

SU  +  2.*PI2«'RaOU»CU-PI2*PI2«RAD*SU-PI2*PI2*DEL2*SUV 

CU-PI2*RADV*SU-PI2«PI2*DEL2*CUV 

SU+PI2*RADV*CU+PI2'!'PI2*DEL2*SUV 

CU-PI2*PI2*DEL1*CV  +  PI2<'PI2«DEL2*CUV 

SU-PI2*PI2*DEL1*SV-PI2*PI2*DEL2*SUV 

>CS1>SC1^SS1/I3S>JBS^K3S*S>0bGSWA^JA,KA) 

^CS2^SC2*SS2>IBU^Jbo>KBu»S^DijGUW3^Jti»KB) 

^CS3fSC3iSS3/IBV,JBVfKBV>S,DBGV,IC>JC^KC) 

*PCS^PSC>RSS/IROiJRO»KRG»S>RO»IR>JR*KR) 

RrS>ROU>RDV»RCC/RCSiRSC^RSS> 

Oi JRO/KRG/ 1R> JR/KR ) 


LUAT (J-l)*RSC( l,J,K) 
LG6T( J-l )*RSS (  1/ J^  K) 
L3iT(J-l)'PCC(I>J*K) 
LGaT(J-1)*RC5(I> J»K) 
RaSU»ROUU^RGUV,TCC>TCS*TSC*TSS* 
0, JkG,KRG> IR/ JR*  KR ) 
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DO  2  J»1*JR0 
DO    2    K«1#KR0 

TCC(I#J»K)«FLaAT(I)»kCC(I+l>J*K) 
TCS(I*J*K)«FLQAT(I)*PCS(H-1*J/K) 
TSC(I*J>K)»FLaAT(I)*R5C(I+l*J>K) 
TSS(I>J*K)«FL0AT(I)*R3S(I+1»J>K) 
CALL  DFSUB(S/ROSS»ROSU»ROSV,TCC>TCS>TSC^TSS> 
H  IROi> JRD* KRUi IR* JR^KR) 

DO  3  I-IWRO 
DO  3  J-1»JR0 
DO    3    K«1>KR0 

TK--FL0AT(K-1 )*FLCaT (K-1) 
TCC{  I> Jf K)«TK*RCC (I» J>K) 
TCS{  1W»K)-TK*RCS  (I>J>K) 
TSC(  1, J*K)«TK*RSC(I* J*K) 
TSS(  I, J^K)«TK*RSS (I> J*K) 

CALL  FSU3(TCC>TCS*TSC>TSS*IR0/JR0*KR0,S>RDVV>IR>JR>KR) 
R0VV»R0VV»PI2*PI2 

COMPUTE  RAX,ZAX,RAX,ZAXV,RAXVV^2AXVV 
KRAXl  «KRAX-1 
KZAXl  «KZAX-1 

CALL  FOURIER(RACfRAS> KRAXl, RAX^CSVZ^SINV) 
CALL  F0URIER(ZAC,ZAS,KZAX1>ZAX,CSV2>SINV) 
DO  ^  K»1>KRAX 
TXC(K)=  FLOAT  (K-l)«(<  AS  (K) 
TXS(K)»-FL0AT(K-1)*PAC(K) 

CALL  F0URIER{TXC,TXS,KRAXl>RAXV,CSV2fSINV) 
RAXV-?I2*RAXV 
DO  5  K»1*KZAX 
TXC(K)»  FL0AT(K-1)*ZAS(K) 
TXS(K)  "-FLOAT (K-1 )*ZAC(K) 

CALL  F0UKIER(TXC,TXS/KZAX1,ZAXV>CSV2>SINV) 
ZAXV«PI2*ZAXV 
DO  6  K=1,KRAX 
TK--FL0AT(K-1  )» FLOAT (K-1) 
TXC(K  )»TK*RAC  (K ) 
TX5 ( K) «TK»RAS ( K) 

CALL  FOURIER  (TXCTXS,  KRAXl, RAXVV,CSV2*S  I NV) 
RAXVV'PI2*PI2*RAX\/V 
DO  7  K«liKZAX 
TK«-FL0AT(K-1 ) ♦FLOAT (K-1) 
TXC (K  )»TK*ZAC (K  ) 
TXS(K)«T.^*ZAS(K) 

CALL  FDJRIER(TXC,TXS,^2Axl>ZAXVV,CS\/2,SIr^;V) 
ZAXVV»PI2*PI2*ZAXVV 
R      «RAX*S*RD» (RGB-RAX) 
RS     • (RO+S«RDS)*(ROB-RAX) 
RU     »S*R0U*(RC6-PAX)+S*RC*R0BU 
RV     •RAXV+S*ROV*( ROB-RAX )+S*RD* (ROB V-RAXV) 
DRSS   «(2.0*R0S*S*R0SS)*(Ra3-RAX) 
DRSU   •(R0U  +  S*RCSU)»(R3  8-RAX)  +  (R0+S*R0S)*'R0BU 
DRSV   ■«(ROV+S*RCSV)*(RJB-RAX)+(RO+S*ROS)*(ROBV-RAXV) 
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DRUS 
DRUU 
DRUV 

DRVS 
DkVU 
DRW 

Z 

zs 

zu 

IM 

DZSS 

DZSU 

DZSV 

DZUS 

DZUU 

DZUV 

DZVS 
DZVU 
OIMM 

RZSU 

DRZSU 

DRZSLi 

DRZSU 

IF(EP 

IF  (EP 

Dl 

DIS 

DIU 

OIV 

DJAC 

DJACS 

DJACU 

DJACV 

D2 

UVUV 

VSJV 

SUUV 

VSVS 

SUV5 

SUSU 

OUVUV 

DUVUV 

DUVUV 

DVSUV 

DVSUV 

DVSUV 

DSUUV 

DSUUV 


S  = 

u- 
v  = 

.L 

.G 


DRSU 

S*ROUU 

S*RQUV 

S»RO*R 

DRSV 

DPUV 

RAXVV+ 

S  ♦  R  Q  ♦  ( 

ZAX+S* 

(RG+S* 
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